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Gliding descent in autorotation is a maneuver used by helicopter pilots in case 
of engine failure. It requires considerable skill, and since it is seldom practiced, it 
is considered quite dangerous. In fact, during certification, a region of low altitude 
and low velocity (the H*V restriction zone) is established where it is considered 
impossible to make a safe descent. 

The landing of a helicopter in autorotation is formulated as a nonlinear op- 
timal control problem. A simplified point-mass model of an OH-58A helicopter 
is used in the study. The model considers as its states the helicopter vertical and 
horizontal velocities, vertical and horizontal displacements (from the point at which 
engine-failure occurred), and the rotor angular speed. It provides an empirical ap- 
proximation for the induced velocity in the vortex-ring state. The cost function of 
the optimal control problem is a weighted sum of the squared horizontal and vertical 
components of the helicopter velocity at touchdown. The control (horizontal and 
vertical components of the thrust coefficient) required to minimize the cost function 
is obtained using nonlinear optimal control theory. 

A unique feature in the present problem formulation is the addition of path 
inequality constraints on both the control and the state vectors. The control variable 
inequality constraint is a reflection of the limitation on the rotor thrust coefficient. 
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The state variable inequality constraint is an upper bound on the vertical sink-rate 
of the helicopter during descent. “Slack” variables are employed to convert these 
path inequality constraints into path equality constraints. The resultant two-point 
boundary-value problem with path equality constraints is successfully solved using 
the Sequential Gradient Restoration Technique. 

Optimal trajectories are calculated for entry conditions well within the H-V 
restriction curve, with the helicopter initially in hover or in forward flight. Solu- 
tions are essentially discontinuous, i.e., they consist of variable subarcs which are 
connected at suitable corners. Subarcs are those which satisfy either the Eulerian 
equations, the upper bound on the thrust coefficient, or the bound on the maximum 
rate of descent. The optimal solutions exhibited similar control techniques as are 
used by helicopter pilots in actual autorotational landings. The results indicate the 
need to drop the collective pitch immediately after engine failure. During the land- 
ing flare phase, the thrust vector is rotated to the rear in order to decelerate the 
forward motion of the vehicle. The stored rotational energy of the rotor is traded 
for additional thrust to cushion the landing at touchdown. The study indicates 
that, subject to pilot acceptability, a substantial reduction could be made in the 
H-V restriction zone using optimal control techniques, thus providing a benchmark 
for comparisons with other control techniques. 

These optimization techniques could also be used to: 

(1) help instruct pilots on good autorotation technique. 

(2) reduce the risk/time/effort involved in establishing the H-V restriction zones 
by flight tests. 


(3) provide an objective comparison of the autorotation capabilities of different 
helicopter models. 

(4) assess the influence of vehicle parameters on autorotation during preliminary 
design. 
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Abstract 


Gliding descent in autorotation is a maneuver used by helicopter pilots in case 
of engine failure. It requires considerable skill, and since it is seldom practiced, it 
is considered quite dangerous. In fact, during certification, a region of low altitude 
and low velocity (the H-V restriction zone) is established where it is considered 
impossible to make a safe descent. 

The landing of a helicopter in autorotation is formulated as a nonlinear op- 
timal control problem. A simplified point-mass model of an OH-58A helicopter 
is used in the study. The model considers as its states the helicopter vertical and 
horizontal velocities, vertical and horizontal displacements (from the point at which 
engine-failure occurred), and the rotor angular speed. It provides an empirical ap- 
proximation for the induced velocity in the vortex-ring state. The cost function of 
the optimal control problem is a weighted sum of the squared horizontal and vertical 
components of the helicopter velocity at touchdown. The control (horizontal and 
vertical components of the thrust coefficient) required to minimize the cost function 
is obtained using nonlinear optimal control theory. 

A unique feature in the present problem formulation is the addition of path 
inequality constraints on both the control and the state vectors. The control variable 
inequality constraint is a reflection of the limitation on the rotor thrust coefficient. 
The state variable inequality constraint is an upper bound on the vertical sink-rate 
of the helicopter during descent. “Slack” variables are employed to convert these 
path inequality constraints into path equality constraints. The resultant two-point 
boundary-value problem with path equality constraints is successfully solved using 
the Sequential Gradient Restoration Technique. 
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Optimal trajectories are calculated for entry conditions well within the H-V 
restriction curve, with the helicopter initially in hover or in forward flight. Solu- 
tions are essentially discontinuous, i.e., they consist of variable subarcs which are 
connected at suitable corners. Subarcs are those which satisfy either the Eulerian 
equations, the upper bound on the thrust coefficient, or the bound on the maximum 
rate of descent. The optimal solutions exhibited similar control techniques as are 
used by helicopter pilots in actual autorotational landings. The results indicate the 
need to drop the collective pitch immediately after engine failure. During the land- 
ing flare phase, the thrust vector is rotated to the rear in order to decelerate the 
forward motion of the vehicle. The stored rotational energy of the rotor is traded 
for additional thrust to cushion the landing at touchdown. The study indicates 
that, subject to pilot acceptability, a substantial reduction could be made in the 
H-V restriction zone using optimal control techniques, thus providing a benchmark 
for comparisons with other control techniques. 

These optimization techniques could also be used to: 

(1) Ms> instruct pilots on good autorotation technique. 

(2) retftice the risk /time/effort involved in establishing the H-V restriction zones 

by flight tests. 

(3) provide an objective comparison of the autorotation capabilities of different 
helicopter models. 

(4) assess the influence of vehicle parameters on autorotation during preliminary 
design. 
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Chapter 1 
Introduction 


§1.1 Background and Motivations of Research 

The unique autorotation capability of the helicopter is an inherent safety feature 
which is heavily relied upon during power failure emergencies. However, the au- 
torotation maneuver, which places great demands on pilot skill, is an unfamiliar 
task to. most helicopter pilots. The vulnerability to power failure has received re- 
ne'wek mterest tately, due to the increased military tactical emphasis on nap-of-earth 
(NOE) operations which require helicopter flight within the restricted areas of the 
height-velocity curve. 


Fig. (1.1.1) shows a typical height-velocity restriction diagram. The “avoid" area 
of this curve defines a region of height and speed combinations from which a given 
helicopter, operated by a pilot of “average" skill, cannot make a safe landing should 
the power source fail. Outside of the avoid area, the pilot has sufficient leeway to 


bldUe 


height 


dllU All iv mamvom 


snsrgy in ths rotor to arrest the 


sink rate at touch down, and thus accomplish a safe landing. 
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While the frequency of emergency autorotative landings has decreased over the past 
several years due to improvements in the reliability and maintenance of helicopter 
engines, the percentage of unsuccessful landings resulting from emergency autoro- 
tations has remained high. United States Army Safety Center accident statistics 
reveal that at least 27 percent of all autorotative landings involving single engine 
helicopters result in some degree of vehicle damage or personnel injury j Reference 
lj. Fig. (1.1.2) shows that most of these emergency landings are related to engine 
failure. 

§1.2 The Autorotation Maneuver and Related Research 

The transient dynamics and control of a helicopter after engine failure have been 
studied both analytically [2,3 . and experimentally 14-10 . The immediate and obvi- 
ous effects of power loss are rotor rpm decay and out-of-trim rotational accelerations 
(notably, left yaw). From the standpoint of minimizing rotor rpm decay, the collec- 
tive pitch must be reduced immediately. This is especially true when collective pitch 
and consequent torque are high. Heavy weight, high altitude, hover and vertical 
ciHmb ate therefore critical conditions. 

Typically, the maneuver of the helicopter, from pilot recognition of engine failure 
to touchdown, can be divided into three phases. The entry phase consists mainly of 
the arresting angular motion of the vehicle and main-rotor rpm decay. During the 
steady-state descent phase, air flows upward through the rotor disk. The increase 
in angle of attack on the rotor blades offsets the reduction in the collective pitch 
angle. Total aerodynamic force is increased and inclined forward so equilibrium is 
established. Potential energy of the vehicle can also be traded for kinetic energy in 
order to attain desired steady-state descent airspeed that correspond to minimum 
sink rate or minimum descent angle. 
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To perform the final phase of an autorotative landing, the pilot must reduce airspeed 
and sink rate just before touchdown. Both of these actions can be accomplished by 
moving the pilots cyclic control stick to the rear. The rearward oriented rotor disk 
allows a larger volume of air to flow through it, resulting in an increase in the total 
lifting force. The increased aft-directed thrust will reduce both the airspeed and 
sink rate. Kinetic energy of the vehicle has been converted into lift (in the forms of 
profile and induced power losses), as well as rotor energy. Finally, the collective is 
raised to convert the stored rotor energy into lift which further cushions the landing. 

Various methods and devices have been proposed to improve helicopter autorota- 
tional characteristics. One passive autorotation augmentation concept is to store 
energy in the helicopter main rotor by using blades with high inertia. Flight demon- 
stration of the concept, the High Energy' Rotor System (HERS), was conducted by 
Bell Helicopter and documented in references [11,12]. In addition to reducing the 
H-Y restriction curves, the HERS can also provide increased maneuverability and 
performance. However, the high payload weight penalty makes HERS unattractive 
for all single engine helicopters. 

Active autorotation augmentation concepts have also been explored. References 
[1,13 and 14] list tip jets, flywheel and auxiliary turbines as the three most promising 
concepts that can provide an additional source of energy to the system with payload 
weight penalty of only 3 to 8 percent. Based upon simulation results, these authors 
conclude that the autorotative characteristics of single engine helicopter can be 
substantially improved. 

In comparison with the concepts of active energy addition and passive energy stor- 
age, the concept of optimal control management as a means of improving autoro- 
tation characteristics of a helicopter has received relatively little attention. Here, 
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improved autorotation performance is achieved only by the management of available 
energy. No supplemental energy is used. 

Johnson '15’ used nonlinear optimal control theory to study the autorotative descent 
and landing of a helicopter in hover. He found that the optimal descent is purely 
vertical. A comparison of the optimal control procedure with flight tests showed 
sufficient correlation to verify the basic features of the mathematical model used. 

§1.3 Objectives of Research 

The primary objective of the research is to study the possible reduction in height- 
velocity restrictions for the autorotational landing of a helicopter using optimal 
control techniques. 

A secondary objective is to develop numerical optimization algorithms that incor- 
porate the practical constraints that are involved in executing the maneuver. 

§1.4 Outline of Thesis 

In Chapter 2, the control of a helicopter after engine failure is formulated as an 
optimal control problem using a simplified point mass model representing an OH- 
58A helicopter equipped with high inertia blades. The formulation contains path 
inequality constraints, reflecting limitations on the rotor thrust coefficient and ver- 
tical sink rates acceptable to pilots. 

In Chapter 3 the numerical optimization algorithms used to solve the problems 
posed in Chapter 2 are described. 
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In Chapter 4, results obtained for optimal autorotative landings of the modeled 
helicopter with entry conditions both inside and outside of the height-velocity re- 
striction curve are presented. These results are compared with those obtained from 
the HERS flight tests [12 , that used a similar helicopter to the one modeled for the 
research. 

Finally, in Chapter 5 w r e discuss the potential usefulness of the optimal procedure 
in the reduction, or even elimination, of the height-velocity restriction curve. Areas 
of further research are also recommended. 

The major contribution of this research is in the formulation of a general optimal 
autorotative descent problem with path inequality constraints, reflecting limitations 
on the thrust coefficient and vertical sink rate. This problem was successfully solved 
using the Sequential Gradient Restoration (SGR) technique. 

In the course of the research, the potential usefulness of two other numerical op- 
timization techniques was identified and algorithms developed. The Combined 
Parameter and Function (CPF) optimization algorithm extends the capability of 
FCNOPT [16' to include an unknown parameter vector in the formulation of the 
optimal control problem. The algorithm SECOND computes neighboring feedback 
control law's for optimal control problems with path equality constraints. 

In addition to its potential usefulness in the reduction of height- velocity restriction 
for helicopter flight, the optimal control procedure can also be used for: 

(1) assessing the influence of basic parameters on the helicopter autorotation 
characteristics during the preliminary design process; 

(2) reducing the time and risk involved in the establishment of the H-V restric- 
tion curves during the helicopter certification process; 
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and (3) providing an objective comparison of the autorotation capability of 
different helicopter designs. 


Chapter 2 

Problem 

Formulation 


In this chapter, the landing of a helicopter after engine failure is formulated as a 
dynamic optimization problem. The assumptions made in the derivation of the 
dynamic model are stated first. The cost function of the optimal control problem 
is formulated as a weighted sum of the square of sink rate and forward speed at 
touchdown. Path inequality constraints, reflecting limitations on the thrust coeffi- 
cient and sink rate, are then introduced. Finally, the end conditions are added to 



§2.1 The Need to Simplify 

The solution of a high-order nonlinear optimal control problem is a difficult task. 
Practical engineering problems, like this one, need to be simplified before current 
optimization techniques can be applied. Practical considerations suggest the use of 
an approximate mathematical model of low order which can describe the dynamic 
system within some tolerable degree of error. Solutions obtained from a simplified 
model of the system often provide a good physical understanding of the problem. 
More accurate models can then be used, if necessary, to include secondary effects 
which were ignored in the simplified model. 
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§2.2 Basic Assumptions 

We simplify the problem here by : 

(A) considering only motion in a vertical plane: 

(B) using a point mass model: 

(C) using a simplified induced velocity model where: 

(1) dynamics of induced velocity are neglected: 

(2) triangular induced velocity distribution is assumed over the rotor disk: 

(3) an empirical determination of induced velocity in the vortex ring 

state is used; 

(D) modeling power losses as follows: 

(1) compressibility and tail rotor power losses have been neglected; 

(2) parasite drag of the fuselage is modeled as an equivalent flat plate area; 

(3) mean profile drag coefficient is assumed constant and independent of the 

angle of attack on the rotor's blades. 

(4) ground effect is neglected; 

(E) neglecting winds and variations in air density. 

Motion in a vertical plane was assumed to keep the number of state variable low 
for the optimization codes. Extension of the point-mass model to three-dimensional 
motions would be straight forward and would include two additional states (heading 
angle and lateral position), and two additional controls (lateral component of thrust 
coefficient and yawrate). 
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Justification of assumption (B) is made through a comparison of the experimentally 
determined steady state sink rate of an OH-58A helicopter in autorotation '12' with 
that computed using the point mass model (cf. Appendix A). Figure (2.2.1) shows 
that the computed steady state sink rate falls between the upper and lower bounds 
of the experimentally determined data. Therefore the neglected pitching motion 
of the helicopter in the point mass model does not enter the dynamic performance 
equations in a significant way. 

The modeling of induced velocity for a helicopter operating in the vortex ring state 
is a difficult task. The approximate formula given by Johnson [15], based upon 
experimental results obtained by Washizu et al [17], is used in the present research 
work. Since the vortex ring state is a condition with high induced power loss, it is 
avoided during autorotation in any event. The error introduced by the approxima- 
tion should be minimal. 

No attempt has been made to incorporate an equation to describe the time rate 
of change of induced velocity. A simple first order inflow lag w'as developed in 
references 19,20 and could be used to refine the present formulation. 

It is well known that the power required to hover near the ground is less than that 
required for hover out of ground effect [21]. However, this performance benefit of 
operating near the ground diminishes rapidly as the airspeed increases [22] . Ground- 
based simulator experiments on the control of a helicopter after engine failure and 
autorotation landings have shown only a minor role played by ground effect in the 
overall autorotational performance of a helicopter [23]. Ground effect is therefore 
neglected in the present study. 

Finally, it is difficult to include the effect of atmospheric disturbances in an open- 
loop optimal control problem. However, neighboring feedback control law's could 
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Figure 2.2.1 A comparison of computed sink rate with 
experimental results of Reference (18) 
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be computed along the nominal optimal path, to convert the open-loop solution 
into closed-loop feedback laws. Deterministic effects of steady wind could easily be 
incorporated in the model. 


§2.3 Equations of Motion 

2.3.1 General Considerations and Coordinate System Used. 

The problem considered here is to find the controls after engine failure to arrive at 
the ground with acceptably small forward and vertical velocity. The helicopter is 
assumed to be in equilibrium level flight at the time of engine failure, with rotor 
speed H, forward speed u, and height ho- 


It is convenient to define the aircraft position from the point of engine failure by 
the coordinates h and x in the vertical and horizontal direction respectively. Fig. 
(2.3.2) shows the coordinate system used. The point at which the engine fails is 
therefore h= 0, and h= ho is the ground. 

Various choices of control variables are possible. One choice is the rotor thrust 
coefficient Ct, and the angle the thrust vector makes with the vertical a. Since a 
is not well defined when Cj becomes very small, and in anticipation of the small 
value of Ct when the collective pitch is lowered after engine failure, it is preferable 
to express the problem in terms of the vertical and horizontal components of Cj: 


Cf, = Ct cos d, 
Ct, = Ct sin a. 


( 1 ) 


The collective pitch control required to obtain this thrust may then be obtained 
from blade element theory as in [24;: 
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DEFINITIONS 

0 Point at vhich engine fails 

T Rotor thrust (lb.vt.) 

D Parasite drag (lb. vt.) 

W Weight of helicopter (lb. vt.) 

x Horizontal distance from point of engine failure 

(ft.) 

h Vertical distance from point of engine failure 

(ft.) 

V Velocity of helicopter (ft/sec.) 

x Horizontal component of V (ft/sec.) 

h Vertical component of V (ft/sec.) 

6 Angle vhich V makes vith the horizon (rad.) 

a Angle betveen thrust vector and vertical (rad.) 


Figure 2.3.2 Coordinate system used 
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h . _ ->».») ... 

(i -!*.<) ’ (2) 

where 0 75 is the rotor collective pitch angle at 75 percent span, while 0 and a are 
the rotor solidity ratio and rotor blade two dimensional lift curve slope respectively. 
The quantities n and A are respectively the advance and inflow ratios defined in the 
tip path plane. With reference to Fig. (2.3.2). the advance ratio fj and inflow ratio 
A are defined as follows: 


V = 
A = 


v cos a -+■ u- sin a 

f Tr ’ 

ti sin a — w cos a + u 


VLR 


(2a) 


Here u' is the vertical velocity, and u is the forward velocity of the helicopter with 
respect to the inertial frame, fl is the rotor angular velocity, and u is the induced 
velocity of the rotor disk. Note that A is defined positive in the positive direction 
of|f|' while p is defined positive in the negative direction of x. 

It is not possible to obtain the longitudinal cyclic pitch control from Cj and a 
without considering the helicopter pitch attitude, and the rotor flapping also. But 
the sign and magnitude of Cj x provides information about the orientation of the 
rotor disk in space. 


2.3.2 Dynamic Equations. 

With reference to Fig. (2.3.2) , vertical and horizontal force balances give: 


m w = mg - T cos a + D sin 6 , 
mil = Tima - D cos 6 , 


( 3 ) 
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Here T is the rotor thrust. The helicopter parasite drag D is defined by an equivalent 
flat plate area f e as: 


= i/>(u 2 T n- ! )/ t . (4) 


The angle 6 which the resultant velocity vector V makes with the horizontal can be 
eliminated from equation (2) by the relationships: 


sin# = 
cos 6 = 


xv 

s/u 2 -r tr 2 ’ 

u 

\/u 2 4- tu 2 


( 5 ) 


Note that 6 is undefined when both v and w approach zero. However, the corre- 
sponding components of parasite drag in the vertical and horizontal direction also 
approach zero under these conditions: 


D sin 6 = - p f e w y/u 2 -r- tr 2 — ♦ 0 

i . — («) 

D cos 6 = -p f t u Vt/ 2 — w 2 — ♦ 0 

2.3.3 The Energy Model. 

A unique characteristic of the helicopter, as opposed to a fixed wing aircraft, is in 
its ability to store energy in the main rotor. The main rotor will accelerate when 
the torque supplied by the engine exceeds the torque required on the main rotor 
shaft. The torque balance equation can be expressed simply as: 

/*ri = -Q, 

= -\p{*R 2 ){nR) 2 R}CQ. 


(ба) 

( бб ) 
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and the energy balance equation of the rotor system is given by: 

i R nn = p s -p R ( 7 ) 

Here I R is the total rotational inertia of the rotor system and Cq is the torque 
coefficient. It can be shown that the torque coefficient Cq is the same as the power 
coefficient Cp [24'. Therefore, we can substitute Cp into equation (6b) for the 
torque balance equation. P§ is the power supplied by the engine and available on 
the main rotor shaft, after losses associated with driving the tail rotor, gearboxes, 
etc. have been deducted. In the event of complete engine failure, power supplied 
to the main rotor is reduced to zero (in fact, shaft power may even be negative 
due to mechanical losses or residual tail-rotor profile losses, etc.). P R is the power 
required on the main rotor shaft to generate lift and propulsive thrust, and also to 
overcome blade profile drag. The induced power (associated with the generation 
of thrust) and propulsive power (power required to overcome parasite drag on the 
fuselage and to accelerate the helicopter forward) can be computed using momentum 
tfep^^^fB&t'drag on the rotor blades must be obtained by blade element theory. 

2.3.4 Momentum Theory. 

In the momentum theory approximation, the rotor affects only the air passing 
through the rotor disc. As the air flows through the rotor disc it experiences a 
velocity increase V perpendicular to the disc (the induced velocity). The thrust 
generated by the rotor is equal to the rate at which momentum is imparted to the 
flow: 

f = -p nR 2 | V - v | 2u ( 8 ) 


since the total velocity imparted to the air flowing through the disc is 2 u (cf. [25’). 
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The momentum power, which is the sum of the induced and propulsive powers is 
simply the rate at which energy is transmitted to the air due to the helicopter flight. 
It is the scalar product of the thrust and the resultant velocity of the flow through 
the actuator disk: 

Pm = T-{V-u) (9) 

The first term in equation (9). Pp — f • V, represents propulsive power. This is the 
power required to accelerate or to climb against parasite drag. The term is negative 
in autorotative flight. The rotor is then like a windmill, extracting energy from the 
air as it sweeps through it. 

The second term, P t = —T • P, represents the induced power required to produce 
thrust. It is always positive since the induced velocity vector is always oriented in 
a direction opposite to the thrust generated. 

Momentum theory cannot account for induced power inefficiencies such as tip loss 
(similar to a rotor with reduced blade size) and those due to non-uniform inflow 
distribution. A tip loss factor of 0.97 has generally been assumed in helicopter 
research and has been neglected in the present work for simplicity. 

For a given thrust, a uniform inflow distribution minimizes the induced power loss. 
A non-uniform inflow distribution raises the induced power by a factor of K tnd and 
the actual induced power requirement becomes: 

P t = -K tnd TP (10) 

where K tnd is the ratio of non-uniform inflow to uniform inflow induced power 
requirements. For a triangular downwash distribution, K ind is given by |( v / 2) 3 , or 
approximately 1.13. 
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2.3.5 Modeling the Induced Velocity. 

The induced velocity v is approximated by Johnson [15' as: 

tv~u = K tnd u h fj f G . (11) 

where the symbols are defined below. 

The ground effect factor fc is taken to be unity in the present study, r is a time 
constant, approximated by [19': 

°.21 

T ~ | a | n 0 

flo is the nominal angular speed of the rotor and is on the order of 37.0 rad/sec for 
OH-58A helicopters. X is the inflow ratio which is of the order of 0.04 (for example, 
^fcot er = 0.039). The value of r calculated using equation (12) is on the order of 0.14 
seconds and may be neglected in our analysis. 

Vh is a reference velocity defined by: 




2pirR 2 

= (13) 


Finally, the induced velocity parameter //is defined as the ratio of the actual 
induced velocity to the reference velocity defined in equation (13). For the deter- 
mination of //, the following expression is used [15]: 
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where the parameters T] and 12 are defined as follow: 


= 


*2 = 


u sin a — w cos a 
Vh 

u sin a — u' cos a 

i?n v 

u cos q — U' sin a 
I'h 

u cos a — w sin a 

i?n v ^ 


(15) 


(16) 


The first expression for //is the familiar momentum theory result. The second 
expression is an empirical approximation for the vortex-ring state (where the mo- 
mentum theory breaks down). The region of roughness in the vortex-ring state is 
defined approximately by [2i\ — 3) 2 + x\ < 1.0 17 . An approximate three dimen- 
sional picture of the variation of // with i\ and X 2 is given in Fig. (2.3.3). 


2.3.6 Profile Power and the Blade Element Theory. 

Accurate descriptions of the profile power require extensive wind tunnel tests to 
determine the effects of thrust coefficient, advance ratio, and angle of attack of 
the rotor's blades on rotor performance. However, in hover and level unaccelerated 
flight, a limited power series expansion of the profile drag coefficient in terms of mean 
blade lift coefficient and advance ratio offers a convenient although approximate 
description of the profile power requirement. 

Following simple blade element theory, the profile power is traditionally referred to 
by an equivalent profile drag coefficient: 


P pro — Cpro P (fli?) 


(17) 



Figure 2.3.3 Variation of induced velocity parameter// 
with i\ and *2 
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The factor c pr0 is the equivalent profile drag coefficient which may be approximated 
by \2 4': 

Cpro = \oc d (1 (^) 2 )(1 - 4.6// 2 ) (18) 

where c d is the mean profile coefficient of the rotor's blades. With the assumed use 
of the NACA 0012 Airfoil on the main rotor’s blades of an OH-58A helicopter, the 
value of Cd may be taken as 0.0087 24.. With values of n and ^ of the order of 
0.15 and 0.06 respectively, both squared terms in equation (18) have been neglected. 
This is acceptable, as the profile power usually represents a small part of the total 
power requirement for helicopters. 

2.3.7 The Energy Equation. 

The total power required on the main rotor shaft is obtained by adding the mo- 
mentum and profile powers together: 

Pr = Pm ~ P i pro 

= TV - K ind T -t7+ P pro (19) 

We next consider the force balance equation of the helicopter in accelerated forward 
flight (see Fig. (2.3.4)): 



mV = T -f mg + D 

(20) 

therefore 

mV • V = f • V t mg- V -r D • V 

(21) 


mVV = ~ (-mV 2 ) 
di v 2 ' 

mg -V = ( mgH ) 

di 


note that: 


and 


( 22 ) 

(23) 
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where H is the height of the helicopter measured in a direction opposite to that of 
g from an arbitrary datum. 

By suitably combining equations (7). (19). (21). (22) and (23) and introducing the 
fuselage parasite power as : 

Ppara — ~ D • 1 

= \^' 3 I. (24) 

(where f e is defined in equation (4) ), the time rate of change of the total energy is 
obtained : 

J t ( mgH -r - ml’ 2 - ±I R D 2 ) = P s - (P, + Ppro + Ppara ) . (25) 

and the corresponding torque balance equation is : 

i R n = -Q, 

= -\p{7r R 2 )(QR) 2 R C P . (25 a) 

where Cp is given by: 

Cp = -‘■Cj'l, (256) 

This energy conservation equation corresponds to the principle that any excess 
power supplied by the engines that is not dissipated by the helicopter is stored as 
internal potential, kinetic or rotational energy. Obviously, the internal energy level 
of the helicopter can only increase if the engine power supplied on the main rotor 
shaft exceeds the total power required. This excess power may be used to climb, to 
accelerate, or to increase rotor speed. 

Conversely, in the event of engine failure, the total power or energy will decrease 
at a rate which depends on the helicopter’s airspeed, main rotor thrust, angle of 
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attack and rotor RPM. The pilot's task during autorotative flight is mainly one 
of Energy Management. This can be achieved through control of the thrust 
vector during descent with the desired result that the aircraft can be landed at a 
desired (achievable) location with small vertical and forward speeds. At the same 
time, during the deceleration phase, the pilot must prevent the main rotor from 
overspeeding which would lead to unacceptable blade centrifugal stresses. This 
task is by no means easy. 

2.3.8 Kinematical Relations. 

The kinematical relations needed in the formulation of the optimal control problem 

are given simply by : 


h = w 

(26) 

x — u 

(27) 


Note that these relations are coupled only one way to the dynamical relations. Since 
h — hQ is a hard terminal constraint on the optimal control problem, equation (26) 
is always needed in its formulation. Equation (27) may however be removed, unless 
there is also a hard constraint on the terminal horizontal distance (as in the case 
where the helicopter is forced to land at a particular spot, perhaps due to terrain 
considerations). The removal of equation (27) will reduce the order of the problem 
and will facilitate the numerical solution. Information on the horizontal distance 
travelled may be obtained through the forward integration of equation (27) after 
the optimal time history of the forward speed has been found. 
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The efficiency and rate of convergence of numerical optimization methods depends 
critically on the scales used for the variables involved. This is especially true in 
nonlinear problems '26 . A “well scaled" problem is one in which similar changes 
in the variables lead to similar changes in the cost function. Now consider a typical 
situation where the engine of the helicopter fails while it is cruising at a forward 
speed of 40 Knots and at an altitude of 400 ft . The magnitude of the thrust 
coefficient Cj used just before engine failure is of the order of 0.003. Rotor speed 
before engine failure is 354 rpm . The different units used by the state/control 
variables, and the range of magnitude that these variables will assume in subsequent 
autorotative descent flight, clearly indicates the need to normalize and to scale. 


The equations of motion may be non-dimensionalized using the quantities flo and 
R. Here flo is the nominal angular speed of the rotor before engine failure and R 
is the radius of the helicopter's rotor. Scaling factors of 10, 100 etc. are used for 
convenience. Non-dimensionalized and scaled quantities for the time, states, and 
controls used in the analysis are defined as follow: 


(a) Normalized time: 




flo f , 
100 1 


(28) 


Hence, one unit of r corresponds to about 16 rotations of the rotor. 


From here onward, the notation ( )' will be used to represent differentiation with 
respect to the normalized time r, where: 



d_ 

dr 


100 d 
flo dt 


(28a) 


(b) Normalized states: 

- ( w \ 
11 _ 'o.olnoJ^ , ' 


( 29 ) 
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and 


*2 = {; 


o.oin 0 i? 

_,n 

X3_( no ) ’ 

h 

X4 UO R*' 

Xh = (lOR^ 


h 


(c) Normalized controls: 


u! = 10 3 C Tt , 

U2 = 10 3 C7' 1 , 

therefore : 10 3 Cr = (tij 2 + t* 2 2 ) * . 


(30) 

(31) 

(32) 

(33) 


(34) 

(35) 

(36) 


The effects of these normalizations are first to convert the time, state and control 
variables into dimensionless quantities, and second to scale them so that they all 
have order of magnitude one. 


If in addition, we also define the following dimensionless constants for the system: 

10 4 ff 


go = 
m 0 = 
/ = 
»o = 


n 0 2 R 

10 pnR* 
m 

U 

20? ri?2 
pnR* 

10 1r 


(37) 


Co = |frc rf (10 3 ) 

I Kind 

y/2000 

po = 0.01(^2000) 


The resultant dimensionless equations of motion are then given by: 
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(a) Dynamical relations (see equations (3) and (7)): 

xi = go ~ rno(u)X3 2 fx \\X\ 2 - xi 2 ), (39) 

x-i = mo(u2^3 2 - fXi^X] 2 — X 2 2 )- (40) 

X 3 ' = -*0J3 2 (C0 •+■ A\/ui 2 T t<2 2 ). (40) 

(b) Kinematical relations (see equations (26) and (27)): 

X 4 — O.lxj . (41) 

x 5 ' = 0.1X2. (42) 


(c) Supporting expressions: 

(cl) the inflow and advance ratios are (see equations (2a)): 


x u sin a — w cos a u 

= M + mi' 


= 0.01 


X 2 U 2 - X\U\ 

X^yUl 2 ■+ t/2 2 


u cos q + tv sm a 
M ™ 1 


+ kofl{ui 2 "t- U2 2 ) 4 . 


= 0.01 


r n 

X2«l + X]t*2 
X 3 \/u ] 2 + U 2 2 


(43) 


(44) 


(c2) the induced velocity parameter // (see equations (14)-(16)): 


X2U 2 -X\Ui 

x\=Po s , 

(45) 

X3(t/j 2 + «2 2 ) 4 


X 2 U 1 -f X\U2 
x 2 = P0 s • 

(46) 

X3(tli 2 -+ U2 2 ) 4 


where once again, the value of the induced velocity parameter //, is determined 
from the expressions given in equation (14): 


2.6 Terminal Constraints and Initial Conditions 


29 


Jl =\ 1 -°/\ // (*2 + (*1 + //)*)* if ( 2*1 ^ 3) 2 + *1 > 1.0; 

( X](0.373xj — 0.598x2 — 1.991). otherwise. 


§2.5 Cost Function 

The optimization problem is to arrive at the ground with small vertical and horizon- 
tal velocities subject to maintaining acceptable conditions during tha autorotative 
descent. The cost function, or performance criterion of the problem can therefore 
be taken as the weighted sum of the squared normalized sink rate and forward speed 
at the time of touch down: 

l,’ + W.V) (47) 

Here H', is the ■weighting function of normalized horizontal speed relative to vertical 
sink rate. Acceptable vertical sink rate at touch down that is compatible with the 
shock absorption capability of typical landing gear design is of the order of 8 fps 
jllj. A reasonable value for forward speed at touch down is 3 knots (this is the 

» v • *■ •••>•• x/... -• _ ..... .... -....y 

average horizontal airspeed at touch down for a series of autorotative descent tests 
on the HERS helicopter {'ll]) . A suitable value of W x is therefore: 

W z = ( ) 2 , 

K 3 x 1 . 688 ; ’ 

= 2.5 . (48) 


§2.6 Terminal Constraints and Initial Conditions 

The helicopter is assumed to be in equilibrium level flight at the time of engine 
failure, with rotor speed fio, forward speed u o, and height of ho. The position of 
the helicopter after engine failure is defined wdth respect to the point at which the 
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engine failure occurred. The coordinate system used is defined in Fig. (2.3.2) where 
h is measured in the downward direction. Therefore, the initial conditions of the 
state variables are: 

£jo - 0, 

£20 = “Or 

£30 = (49) 

£40 = Or 
and £50 = 0 

where uq is defined to be ( 0 o“&n 0 ) • 


The terminal constraints of the optimization problem include: 

£ 4 / = hj , 

*» 1 - d i- 


(50) 


where: 



(51) 


Note that, while the first equation of (50) is always needed, the second equation is 
used only when there is a hard constraint on the terminal horizontal distance (to 
land at a horizontal distance of df ft away from the point at which engine failure 
occurred). 


§2.7 Path Constraints 

The equivalent profile drag coefficient of the rotor (cpr 0 , as defined in equation (18) 
) increases sharply when the thrust coefficient exceeds the rotor stall limit (£z . ) ltaU . 
The immediate effect of this increase in profile drag is a drop in the rotor speed. This 
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drop in rpm causes an increase in th«- angle of attack on the rotor’s blades and will 
ultimately lead to rotor stall and the instability associated with it. This limitation 
on the thrust coefficient requires a path inequality constraint in the optimal control 
problem. 

A typical value of ) t t a ii for the OH-58A helicopter is 0.15. This value is used 
in the present study. The path inequality constraint, and its non-dimensionalized 
form are: 

or 7.2 > v/«i 2 + u 2 2 (52) 

where: 

(10 Cx)itall — I® ( 

o 

= 10 3 x 0. 15 x 0. 048, (53) 

= 7.2. 

This inequality constraint can be converted to a path equality constraint as shown 
below: 

Since (C j ,) 2 > (uj 2 + U2 2 ) 

where Cj, is equal to 7.2. 

therefore Cj t - (uj 2 + U2 2 ) - t/3 2 = 0. (54) 

where U3 is a “slack variable” or artificial control that has been introduced to convert 
the inequality constraint into an equality constraint. 

The upper bound on the vertical sink rate as an additional path inequality constraint 
will be discussed in Section (3.3). 
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§2.8 Further Time Normalization 

The optimal control problem that has been posed thus far is one with an unspecified 
terminal time Tf. The problem may be converted into one with specified terminal 
time through the following (further) normalization of the dimensionless time r: 

f = T - (») 


The transformation (55) converts the independent variable from r to £ where £ now 
varies from 0 to 1. This transformation introduces into the problem an additional 
unknown parameter Tf that has to be optimally selected. We shall from here onward 
denote the differentiation wfith respect to £ by: 


( ) V = — 

1 ' dC 


= r f ( )'• 


(56) 


§2.9 Final Form of Helicopter Optimization Problem 

We are now in a position to write dotvn the final form of the helicopter optimization 
problem, Let: 


X = (xi x 2 x 3 n x 5 ) , (57) 

U = (uj u 2 u 3 ) T , (58) 


7f = {T f ). 


(59) 


denote the state, control, and unknown parameter vectors of the problem. 
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The problem is to find J7(£) and rr to minimize: 

/= + = *(X,), 

subject to : 

(1) equations of motion ( X v = f ): 

x i V = T f{9o ~ m 0 (u ix 3 2 - fx i \ xi 2 -r x 2 2 )), 

j 2 V = r / mo(u 2 x 3 2 - /x 2 v'xi 2 + x 2 2 ), 

*3^ = -T/*0 X 3 2 (C0 + A\/ U] 2 + u 2 *), 
x 4 v = O.lr/Xj, 

X 5 r = O.lr,^. 

(2) the initial condition of A' is given by: 

Ao = (0, uo, 1,0, 0) r . 

(3) path equality constraint (S(X,U,H) = 0): 

(Isi+i!!)* . o. 

C T. 

(4) terminal constraints (t^( A/,?r) = 0): 


(60) 


(61) 


(62) 


(63) 


x 4f -h f = 0. 
X s f - df = 0. 


(64) 


In the next chapter, a gradient-type numerical optimization technique that can be 
used to solve this problem is described in detail. It requires the calculation of 
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(5x3 matrix). ^4 (5x5 matrix), (5x1 matrix), ^5 (1x3 matrix), (l>5 
' o A OTi oL o\j 

matrix) and (2x5 matrix). Detailed expressions of these matrices are given in 
Appendix (B). 


Chapter 3 

Numerical 
Optimization Techniques 


In the previous chapter, the landing of a helicopter after engine failure was for- 
mulated as a nonlinear optimal control problem with path equality constraints. 
Numerical optimization algorithms that can be used to solve this problem are de- 
scribed in this chapter. 

, . .. . Vv ; . .•l' ; . . , \ • 

• 2 * : ; • - - . •%, . 

This chapter begins with a review of algorithms for solving optimal programming 
problems with bounded controls and/or states. A combined function and parameter 
optimization algorithm is then described. It is an extension of the ordinary gradient- 
type numerical algorithm (FCNOPT) [16], to handle the presence of an unknown 
parameter vector. In Section (3.3), we describe the Sequential Gradient Restoration 
algorithm which can be used to solve optimization problems with nondifferential 
path equality constraints. Several transformation techniques are then presented 
that convert problems with path inequality constraints to problem with equality 
constraints. The chapter ends with a description of an algorithm that can be used 
to compute neighboring feedback control laws for optimization problems with path 
equality constraints. 
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§3.1 Algorithms for Problems with Bounded Controls and/or States 

One of the earliest attempts at numerical solution of optimal programming prob- 
lems with control or state inequality constraints was made by Bryson et al [27-28. 
see also chapter 3 of 29 . Necessary conditions for extremal solutions to program- 
ming problems with an inequality constraint on a function of the control or state 
variables were given. It was shown that, in general, certain terms must be added 
to the Hamiltonian function during the interval in which the solution curve lies 
on the constraint boundary. Furthermore, for inequality constraint functions not 
explicitly involving the control variable, one or more functions of the state and time 
must satisfy equality constraints at the beginning (the entry corner) of a constraint 
boundary interval. The Lagrange multiplier functions are not uniquely defined on 
state constraints. In the work reported in Reference 28, a modified version of the 
steepest-ascent technique was used in the numerical solution of two atmospheric 
entry trajectory problems. An advantage of the method is that improvements in 
the control program are not required for the period on the constraint boundary, 
making possible a more rapid convergence towards the optimal program. However, 
the method requires prior assumptions concerning the number and location of the 
junction points. 

Inequality constraints on functions of control and/or state variables have also been 
treated by several investigators through the use of integral penalty functions [31-33/ 
One such scheme [31-32] uses an auxiliary state variable which is the integral of a 
quadratic measure of the violation of the inequality constraint, which is brought as 
close to zero as is necessary to provide a satisfactorily small violation of inequality 
constraint. Rate of convergence to a satisfactory solution is usually slow [29]. 

McGill [33] developed a generalized Newton-Raphson algorithm based upon essen- 
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tially the same idea of reducing the constraint problem to an unconstrained one 
by the introduction of an additional state. His approach does not require assump- 
tions as to the number and location of junction points. Computational experience 
with one example problem (a modified version of the classical brachistochrone prob- 
lem) suggests that it may be useful for obtaining solutions to the class of nonlinear 
problems with bounds on the state space. 

In an extension of the work given in references [27-28], Speyer and Bryson [34 1 
gave a new set of necessary conditions for solution of an optimal programming 
problem with a state variable inequality constraint. It was shown that unconstrained 
arcs must satisfy certain “tangency" constraints: namely, these arcs must h av& 
zero values of the state variable constraint function and all of its time derivatives 
that do not involve the control function, at both the entry and exit corners of 
the constrained arc. These conditions are satisfied automatically if the necessary 
conditions of reference [27] are used. However, if one uses the “direct-adjoining’" 
approach of reference [35], explicit use must be made of the tangency constraints 
at both corners. 

In reference |36 ia Mehra et al showed that some of the difficulties associated with 
nonlinear programming problems with state variable inequality constraints and sin- 
gular arcs arise due to the exclusive use of control variables as the independent 
variables in the search procedure. They proposed a conjugate-gradient algorithm 
which based its choice of independent variables on the problem constraints. This 
choice could result in different combinations of state and control variables as inde- 
pendent variables along different parts of the trajectory. Four numerical examples 
were successfully solved using this approach. Two of the solved problems had a state 
variable inequality constraint and the other two had singular arcs. The inequality 
constraint and singular arcs were handled in a regular fashion without explicit use 
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of extra necessary conditions of optimality. This is considered to be a special feature 
and an advantage of the method. 

In reference [37), Maurer et al distinguish between two main cases of optimal pro- 
gramming problems with bounded state variables, depending on whether the control 
variable appears nonlinearly or linearly. The distinction arises naturally from the 
fact that a nonlinear optimal control in the first class must be continuous at the 
junctions between the interior arcs and boundary arcs, whereas any linear optimal 
control in the second class is discontinuous in general. Maurer et al exploited the 
fact that the optimal control in the first case must be continuous and satisfy a 
suitable augmented Two-Point Boundary- Value Problem (TPBVP). This consists 
of the basic two-point boundary- value problem of the Ma xi mum principle, and of 
additional differential equations and boundary conditions constructed in such a 
way that all the necessary conditions are automatically fulfilled by the solution. 
The TPBVP’s encountered were solved using the method of multiple shooting [38'. 
Three numerical examples were solved to illustrate the efficiency of the algorithm. 
Once again, predetermination of the number and sequence of subarcs of the optimal 
solution has to be made prior to the initialization of the algorithm. 

Most of the above mentioned methods suffer from the following disadvantages: 

(1) it is not known a priori, whether there will be more than one joining of 
unconstrained and constrained arcs, and where in time these joining will 
occur; 

(2) the discontinuities in the Lagrange variables at the junctions of the 
unconstrained and constrained arcs are, a priori, unknown, and have 
to be guessed and iterated upon. 
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To overcome these difficulties, Jacobson et al 139] suggested the use of a different 
approach. A device suggested by Valentine was used to transform a control problem 
with a state variable inequality constraint into an unconstrained one of increased 
dimension. With a p th order state variable inequality constraint, it can be shown 
that in the transformed unconstrained problem, the p th time derivative of the 
slack variable becomes the new control variable. One feature of the transformed 
problem is that any guessed nominal control results in a feasible state trajectory, 
i.e., the inequality constraint is automatically satisfied. A second feature is that the 
transformed problem exhibits singular arcs which correspond, in the original state 
constrained problem, to arcs lying along the state constraint boundary. 

However, a major difficulty arises with Jacobson’s approach when the number of 
state bounds (r) does not equal the dimension of the control vector (m). In partic- 
ular, if m < r. then one cannot express u as functions of the r slack variables unless 
some of these variables are dependent upon each other. Thus one cannot use the 
appropriate time derivative of the r slack variables as independent new controls. 
Tw o nu meri cal e xamples without the above difficulty were given in [39] . 

• •• • • c ••• 

Instead of replacing the control u in the original problem formulation by appropriate 
slack variables and its time derivatives, one can enforce the state/control bound 
by the addition of path equality constraints to the original problem. These path 
equality constraints are again obtained by the use of Valentine’s device on path 
inequality constraints. In this way, nonlinear optimal programming problems with 
path inequality constraints can be transformed into ones with nondifferential path 
equality constraints. 

The solution to this class of problems with path equality constraints differs from the 
solution without the constraints. The usual approach of backward integration of 
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the adjoint-equations and forward integration of the equations of motion cannot be 
employed here. This is due to the fact that the computation of the Lagrange mul- 
tiplier p(t) (associated with the path constraints S) over the time interval requires 
the simultaneous solution of both the adjoint-equations and the equations of mo- 
tion. Because of the coupling, the total system must be integrated simultaneously 
in either the forward or backward direction. One way of solving the coupled, non- 
linear, two-point boundary- value problem is by the Method of Particular Solutions 
given in references 40-41 

Miele et al ;42’. at about the same time as Jacobson, developed a Sequential Gra- 
dient Restoration algorithm for the solution of optimal control problems with path 
equality constraints. The algorithm made use of the method of particular solutions 
of 40'- 41;. Both the feasibility as well as the efficiency of the algorithm were il- 
lustrated through the solutions of example problems given in 42 and other related 
papers '43-45'. The approach taken in the present work on the optimal autorotation 
trajectory of a helicopter is to combine these two mathematical tools (the Valen- 
tine’s device and the Method of Particular Solutions for the solution of a nonlinear 
TPBVP) for the solution of the constrained optimization problem posed in Section 
(2.9). The Sequential Gradient Restoration method is described in greater detail in 
Section (3.3). 

§3.2 Combined Function and Parameter Optimisation Algorithm 

Direct analytical solutions of dynamical optimization problems are only possible 
when the system equations, the performance index, and the constraints of the prob- 
lem are very simple. One of the more reliable methods for the numerical solution of 
the dynamical optimization problem is the steepest-descent gradient method [30’. 
The method is characterized by iterative steps for improving estimates of the con- 
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trol history, u(<), so as to come closer to satisfying the optimality conditions and 
the boundary conditions. First-order gradient methods usually show great improve- 
ments in the first few iterations but have poor convergence characteristics as the 
optimal solution is approached. Second-order gradient methods have excellent con- 
vergence characteristics as the optimal solution is approached but may have starting 
difficulties associated with picking a “convex” nominal solution {29;. Variable Step- 
size Gradient methods can improve the convergence characteristics of the first-order 
gradient methods through the optimal selection of step-size at each iteration step. 

Further information on the variable step-size gradient algorithm can be obtained 
from [26’. 


In this section, we wish to extend the gradient algorithm for function optimization 
problems [16 J so that it is capable of solving optimal control problems that also have 
an unknown parameter vector. We assume that there are no inequality constraints, 
that the initial time and state are fixed, and that functions of some of the state 


variables and unknown parameters are specified at a given or an unspecified terminal 


time. Thus we wish to consider the following optimization problem: 

min/ = if) + / L(x,u,n,t)dt, 

■ -'SviV •»* Jo 

~ , X — /(x, U, TT,t ), 

( 1 ) 

x(0) = given. 


and tp(£f, tF) = 0. 


Here x(n x 1), u(m x 1), and n{p x 1) are the state, control and the unknown 
parameter vectors. tp(q x 1) are the terminal constraints of the problem. If the 
end-time t f of the problem is unspecified, the problem may be converted to one 
with a specified end-time through the use of the following transformation: 


I 
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The independent variable is now r which varies from 0 to 1 (therefore the subscript 
“1” in equation (1) denotes terminal condition). The final time tj is included as 
one of the components of 7r. L and / are replaced by L = tj L and / = tjf. 

If we adjoin the system differential equations and the terminal constraints to the 
performance index I with multiplier function A(/)(r? > 1) and multiplier n(q x 1) 


respectively, we get: 


J = (4>-r n T x. 


>,♦/! 

JO 


L -h X J ( / — i)j dr . 


Now consider the variation in J due to variations in both the control and parameter 

vectors: ^ 

6J=(-\ t + * 1 ) 1 6 x 1 + \[ Hxdr+i*,)^ 

,« ° ( 3 ) 

+ / \H u Sit + (H z -f X T )Sx]dr , 

Jo 

w’here for convenience, we have defined scalar functions H (the Hamiltonian func- 


tion) and $ as given below: 


H = L-rX T f, 
$ = d > + . 


Therefore, first-order necessary conditions for the extremal solution are given by 


the following relations: 


y = -H z , 

y T = (*«h, 

H u = 0 , 


and f H„dr + ($ ff )i = 0.- 
J o 


In the special case when the unknown parameter W contains the unspecified ter- 
minal time, the corresponding component in the last equation in (5) becomes the 
transversal ity condition of the classical literature (cf. Appendix D). 
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Equations (l) and (5) contain all the conditions needed in the solutions of x(r ) (n x 
1 ), tT(r) (m x 1 ), A(r) (n x 1), if [p x 1 ) and (I(q x 1 ). 

Various methods have been suggested for the solution of this two-point boundary- 
value problem. The method used here may be termed the Impulse Response method 

[29:. 

The method involves making initial guesses for both the control time history u(r) 
and the unknown parameter if. The dynamic equations are integrated forward us- 
ing the given initial conditions. In general, the terminal conditions are not satisfied. 
To make improvements in the feasibility condition, we consider the following (q->- 1 ) 
impulse response functions H u ^\ where i = 0 , 1 , •••, 9 . Note that repre- 

sents the variation in the cost function due to a unit impulse (Dirac Function) in 6u 
at time r, while holding x( 0 ) and if constant and satisfying the dynamic equations. 
Similarly, Ff u ^‘*(r) where i = l,...,q corresponds to variations in the terminal con- 
straints V due to a unit impulse in S u at time r. These impulse response functions 
are given by the following expression: 

a.W(T) = rl'l£„(r)-;AW: T /„(r), (6) 


where 



if t = 0 ; 
if i ± 0 . 


(6a) 


for i = 0 , 1 ,. ..,q. The influence functions A(')(r) are determined from the backward 
integrations of the following adjoint equations using the time histories of x (r), u (r) 
and the value of if from the forward integrations: 


AM = -rMi,(r) - 

I*" 1 ), = (*)/. 


0 ) 

(8) 
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where 

xJjq = (9) 

In this way, the (q-l) impulse response functions H u ( '\ where i = 0,1,. ..,q, for the 
current iteration, can be computed. 

The quantities '/J H^ 3 \r)dT 4- (t L'j) n }, where j = 0,1.. ...q represent the variations 
in either the performance index or terminal constraints due to a unit changes in 
the elements of n, while holding f(0) and tT constant and satisfying the dynamic 
equation. 


Therefore variations in the cost function and the terminal constraints due to small 
variations in the control histories 6u{r) and the parameters are given by the 
following relations: 


and 


6J = [' Hj°>6u{T)dT + \f Hj°\r)dT 4 {rbo) T \f>*, (10) 

Jo Jo 

Svj = f H u M(r)6udT + [ [' H n M(T)dT - (v,)* «jf , ( 11 ) 

Jo Jo 


where j = l.....q 


If w’e adjoin equation (11) to (10) with q Lagrange multipliers n t , where i = l,...,q, 
we have: 


Sj = b J 4- ^ , 

t=i 

= / [H U W + Yti t H u M}6udT 

Jo «= i 

+ i f H^ 0) dr + {tpo) n + f H„ {,) dr + Y'n t (v t )„}6n. 

Jo ,=i Jo ,=i 


( 12 ) 
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If we wish to obtain the largest change in 6 j, we would calculate the gradients 
and [J 0 > H„W(r)dT + (0,)„] (i = 0,l v ..,q) and then direct SU and Su opposite 
to their respective gradients as follow: 

» = ~K .\ /' H,"Ut + (<!»). (' 

• /o ,=i * 

= + (fc), 1- + £*( f(H^) T dr + (v,)„ r )' , (13) 

• ,0 .=i h 

where if* is a (p v p) diagonal gain matrix that controls the stepsize of SW. 


Similarly, 


6v = -A' B ;(/f tt W) r + ^ Mt (^« (,) ) r j, 


(14) 


«=i 


where K u is an (m x m) diagonal gain matrix that controls the stepsize of Su. 

When we make these choices in Sir and Su, the predicted variation in the augmented 
performance index J = J — is: 

« = -A* r j f\Hj°>) T dr + (0o), r + V„,( ['(H.MfdT - ((*)/)) 


(15) 


Here the weighted norm-function A ’k( ) is defined as: 

Kk, \y, = y t k,y. 


( 15 °) 


Similarly, the substitutions of equations (13) and (14) into (11) give the following 
changes in the terminal constraints: 

H'j = - f\H u M)K u (Hj 0 YdT-Y»i [' (H u M)K u (H u W) t dr 
Jo t=1 Jo 

- [ f 1 HJ j) dr + (t/>,) .](*,)! ('{HrWfdr + (tfo), T l 
Jo Jo 


( 16 ) 
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where j = 

One way to simplify expressions (16) and (17) is to “absorb" 1 K u into and A* 

into To do this, we first note that both K u and K v . being positive definite 

matrices, can be expressed as products of their respective “square-roots": 

K„ = SrSj , 

K u = S u Sj . 


( 17 ) 


Let us define the following quantities: 

H [ u 0 ] (t) = Hj 0] (r)S u , 
tffV) = H u {,) (t)S u . 


(18a) 


where i ~ l....,q. 


Therefore an expression for H u {t) is given by 


H u (t) = Hj 0 '>(T)^Y t n,H u M(T) 
1=1 

H„(r) = HJr)S u = 


( 186 ) 


1=1 


Similarly, i/P and (V( )) ff are defined as 

ff< 0 ) (r) = ff n < 0 >(r)S„, 

(19) 

M* = M n s„, 


where* = and j = 1 
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Substitutions of (17)-(19) into (15)-(16) give 

6J = -N ['{SPfdr + (*)/ + £>,( ['(H^/dr * (&)/)j 

Jo f=1 Jo 


1=1 


H'j = - [' ffl j) (M° } ) T dr - Tn t f HL }] (Hl' ] ) T dT 

Jo l=I Jo 

- \[' til j) dr - (vJJ: ['(H^fdr + (&)/] 

Jo Jo 

- V «,! /’ + (r,), j f\aP?iT + (&). T | ■ 

“T Jo Jo 


(20) 


( 21 ) 


where j = 

Equation (21) can be rewritten in a vector-matrix form after the introduction of 
the following notation: 

/I(flxl) = (/ii....,// 4 ) T , (22) 

T 

and g (g x l) = (gi,...,g 9 ) . 


The j th component of g is given by: 

9, = - j' Hl’HftFYfr - I !jf ' H^ir + ( 0 ,-),l \J\si 0> ) T dr + (*)/] (23) 

where j = 0,1,... ,q 

Furthermore, let us define a (q x q) symmetrical matrix Q whose components are 

Qtj = Qji i 

= - /' alWY* - 1 f hI’Ut + (^)j| [\ttPfdT + (v;,), r ) . (J4) 

Jo Jo Jo 


where (»,j) = 1, •, 9 
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Using equations (22)-(24). equation (21) can be written in a more compact form: 

f>V = g-QP- (25) 

If we choose by = -ey, where 0< c <1. the Lagrange multiplier p can be 
evaluated as: 

P = -Q~ 1 {g -r tv) . (26) 


If Q~ l exists (the controllability condition, see also Chapter 2 of 29 ), we can sub- 
stitute the value of p found in equation (26) into both equations (13) and (14) and 
express both bH and bv only as functions of the impulse response functions H u ^(t) 
and H r ^(r) (that have been computed and stored during backward integrations): 


Sir = -Sy f\si 0) ) T dr - (v 0 )r 7 - ['{HPfdr + (0, 

Jo —r Jo 




SU= -S u [(hI°Y j. 


For the next iteration, the uses of new trial values 

tyv-i(r) = tT v ( ~ ) + buy(r) . 
~ H K “h N ? 


(28) 


(where N is the current iteration number) 

will reduce the value of the performance index and bring the system closer to sat- 
isfying the terminal constraints. 

After many iterations of the above procedure, all the terminal constraints are sat- 
isfied to a prescribed accuracy: 


I vM i/.ir) i< < v , 


(29) 
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for j = where c^. denote? the feasibility condition (e.g. c v = 10 -4 ). 


Using this procedure, a feasible path satisfying the constraints is determined. How- 
ever, the performance index can still be decreased further. The procedure now 
concentrates on minimizing the value of the performance index while keeping 
within the accuracy of (29). The procedure is repeated until 


f N\{Hl°y + Y M&l'V.dT < <U 
h 


(30) 


and 


.v; f\a l°Y<ir + (*£ + £>,< oi) 

Jo Jo 

1—1 '*»»* 


where e u and £* represent the optimality conditions (e.g. £„ = £„ = 10 -8 ). 


A numerical algorithm that implements the above procedures has been written. 
The procedure is basically an extension of the program FCNOPT by Bryson ;16 . 


Example Problems 

Two^exaas^ie applications of the Combined Function and Parameter Optimization 
Algorithm are given here. The first example involves a nonlinear first order system 
with a specified control law. The second example considers the optimal landing of 
a helicopter that is initially in hover. 


(1) Specified Control Law Problem [46] 

Consider the following optimization problem: 

1 t 1 

min I = - / (x 2 -l- u 2 ) dt 
“ 2 Jo 


subject to the scalar nonlinear differential relation 


} 
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The initial condition is 
There is no terminal condition. 

In order to exercise the parameter optimizatiobn capability, the control law is spec- 
ified to have the following form 

u = a j , 

where a is an unknown parameter whose value is to be optimally selected. The 
transformed problem becomes 

• r 1 
mm I = - 

a 2 

with the following equation of motion 

f = -i 2 -fox, 

while the initial condition of x(0) remains unchanged. 

The transformed problem now becomes a problem with an unknown parameter a 
but without any control. The problem can be solved using the combined function 
and parameter algorithm and the optimum value of a found is: 

a = -0.10334. 

The minimum value of the cost function is: 

I = 4.5218. 

The value of 1 is slightly larger than the minimum value when u(t) is open: 

Im in 
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x(0) = 10. 


= 4.5108. 
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Time histories of both the state and control of the system, using either u(t) — a x(f) 
or the open-loop optimal solution are given in Figure (3.2.1). 

(2) Optima] Landing of a Helicopter in Hover 

In reference [15], Johnson formulated an optimal, autorotative landing problem for 
a helicopter initially in hover close to the ground. He showed that the optimal 
descent is purely vertical. Under these conditions, the general landing problem 
posed in Chapter 2 can be simplified. 


Imposing a purely vertical descent, the horizontal distance component x, velocity 
component u, and control component Cj x are omitted. The remaining states of the 
problem are the vertical height h. vertical sink rate V and the angular speed of the 
rotor fi. In addition, the induced velocity parameter fj in the vortex-ring state 
and the momentum theory states is given by the following one-parameter family of 


equations: 


fl 



\ v + \A + (y ) 2 ’ 

—V’ (0.373V' 2 - 1.991), 



V ( t ) 2 - 1 


for 0 < V < 1 
for 1 < V < 2 
for 2 < V 


(32) 


Note that the first and third expressions of equation (32) are solutions of the mo- 
mentum quartic, while the second is an empirical expression given by Johnson [15] 
for the vortex-ring state. The dimensionless parameter V is defined as the ratio of 
the vertical sink rate and the induced velocity in hover »/*: 


V = 




( 33 ) 


Figure (3.2.2) shows the variation of // with the parameter V. 


SlaLe 
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Figure 3.2.1 Example Optimal Control Problem 






Figure 3.2.2 Variation of Induced Velocity Parameter // with V 
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If we go through the normalization and scaling processes of Section (2.4), the non- 
dimensionalized and scaled time, states and control of the problem are given by the 
following relations: 


_ _ ,n°/ 

7 ^ 100 ^ 


= (; 


v 


1 0.Olfloi? 

( n \ 

X2 = { n 0 ] ' 

t h \ 

X3 = i R h 


)• 


(34) 


U] = 10 3 Cj z . 

In terms of these normalized quantities, the optimal autorotation from hover is 
formulated as the minimization of the cost function /, where 




The equations of motion are given by 

i\ - go ~ moilui - f 0 x ] , 

7* i 

z"2 = to^2 u i ( 0.01 -r- A'o\ / u if]) , (35) 

I 2 

13 = Xi . 

Here ( )' denotes time differentiation with respect to f and the dimensionless quan- 
tities go , mo, fo ■ co, and ko are as defined in Section (2.4) (see equation (37)). The 
induced velocity parameter /; is a function of V, which is now given by the following 
normalized expression: 

V = po— %=, (35a) 

x 2 V U 1 

where po is also defined in Section (2.4). 

The optimization problem is constrained at the end-time by the “stopping” condi- 
tion: 


X3f -hf = 0 , 


(36) 
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and the following one-sided path inequality constraint: 

< &T.- ( 37 ) 


Further Time Normalization 


Since the end-time of the problem is unspecified, further normalization of the inde- 
pendent variable f is needed. The problem is converted to one with fixed end-time 
by the following change of independent variable: 

e = T f ’ ( 38 ) 


The transformation (38) converts the independent variable from f to £ where £ now’ 
varies from 0 to 1. This transformation introduces into the problem an additional 
unknown parameter T f that must be optimally determined. From here on we shall 
denote differentiation with respect to £ by: 



= ff( )'• 


(39) 


Substitution Technique 

The path inequality constraint (37) can be eliminated if we substitute for t/j the 
new control variable U 2 , where U 2 is defined as : 


Uj = Cx, COS 2 U2 . 


(40c) 


Final Form of Optimal Autorotation Problem 
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The state, control, and unknown parameter vectors of the problem are: 

X = (zi X 2 *3)^ (41) 

t/ = (ti 2 ). (42) 

7T = ( t } ). (43) 

denote the state, control, and unknown parameter vectors of the problem. 

The problem is to find U{£) and i r to minimize: 

/=j[n /!’ = (44) 

subject to: 

(1) equations of motion ( — f )'■ 

x \ V = Tf[g 0 - moX 2 2 (Cj, cos 2 u 2 ) - /oX] 2 i , 

x 2 v = -7yt'ox 2 2 co + {Ct, cos 2 u 2 ) (— 0.01 -r fcoy Cr, cosu 2 /;)_ , (45) 

X 3 V = f/X] . 

(2) the initial condition of X is given by: 

X 0 = (0,1, 0) 7 *. (46) 

(3) terminal constraint (rp(Xf,i r ) = 0): 

*3f~hf = 0 (47) 

The optimization problem (41)-(47) is now in a form that can be solved using 
the combined function and parameter gradient algorithm. Several other example 
applications of the algorithm are given in Appendix C. 


3.3 Sequential Gradient Restoration Technique |42 


57 


§3.3 Sequential Gradient Restoration Technique [42] 

This section explains the Sequential Gradient Restoration (SGR) method given in 
[42 . The SGR technique is one way that problems with bounded control may be 
solved. The technique is employed in the solution of the optimal helicopter landing 
problem posed in Section (2.9). 

The problem is similar to the one in Section (3.2) except that here we consider 
more general path constraints. Specifically, we consider the following problem with 
unknown parameter vector t and path equality constraints 5: 

min/ = <f>(xt,r) + / L(x,u, n,r)dr, 

Jo 

x = /(x, u, tt, r), 
x(0) = given, 

5(x, u, 7r, t ) = 0, 
and rp(xf,H) = 0. 

Where 5(r x 1} are the path equality constraints of the problem. The motivations 

;• ■ •< • • •: 'i ■. ■■ :V. • - 

for this formulation" are: 

(1) Many optimization problems arise directly in the form considered here. i.e. 
with path equality constraints. 

(2) Problems involving inequality constraints can be converted to the present 
scheme through suitable transformations. This statement applies, for 
instance, to the following situations: 

(a) Problems with bounded controls. 

(b) Problem with bounded states. 

(c) Problems with bounded time rate of change of state. 


58 


3. Numerical Optimization Techniques 


(d) Problems where some bound is imposed on an arbitrarily prescribed 
function of the parameters, the control and the time rate of 
change of the state. 

The transformation techniques employed are of the Valentine-type [41]. The trans- 
formation is performed by augmenting the dimension of the control space, with 
or without augmentation of the state space. In the process, non-differential path 
equality constraints relating the original control and the slack control are produced 
that must be satisfied by the solution of the original optimization problem. Four 
example problems with bounded control, state or time rate of change of state are 
presented in Appendix (G). 

3.3.1 First Order Conditions 

If we adjoin the system differential equations, the path equality constraints and the 
terminal constraints to the performance index I with multiplier functions A(r)(n x 
1) and p(~)(r > 1), and multiplier p(q x 1) respectively, we obtain the following 
augmented performance index J 

J = (o + li T t'.'jj + f [L + X T (f - i) + p T S'\dr , (2) 

Jo 

= (A T x) 0 - (*- X T x) 1 + f {\ r x + H)d-, (3) 

Jo 

where, for convenience, w r e have defined scalar functions H (the Hamiltonian func- 
tion) and $ as given below: 

j>, r) - L(x,u,n,r) + A r /(£,tT,^r) -r p r S(x,u,7r, r) , (5) 

= <£(*/,*) + At r 0(£ / ,7r) . (6) 
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Now consider the variation in J due to variations in both the control and the 
parameter vectors 

6J = (X T Sx) 0 - (~X T - - ! / 1 H*d? -r ($„),]£* 

Jo 

+ / (tf z + A r )6f dr. 

Jo 

Therefore, first-order necessary conditions of the problem are given by the following 
relations: 

x = f(x,u.x,r ) , 

S(x, u, 7f,r) = 0, 

• T 

X =-H x , 

H u = 0, 

( 4 ) 

x(0) = given, 

i^), = (*»)i, 
tj;{xf,i r) = 0, 
and (#ir)i + f H*dr = 0. 

Jo 

Summarizing, we seek functions x(r), u(r), f and multipliers A(r), p[r) and n which 
satisfy the first order necessary conditions (4). 

3.3.2 Approximate Methods 


In general, the problem posed in Section (3.3.1) is nonlinear and can only be solved 
iteratively. In this connection, we define two scalar functionals 

p = C N\x-f\dr + f\\S}dr + JV|(0),j, (7) 

Jo Jo 


and 


Q = /' N\X+Hj\dr + f N\Hl}ir + N\ (*?),+ /’ Hf*] + A'[(A-*f)i! . 
Jo Jo Jo 


(8) 
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where the norm-function N ] was defined in equation (15a) of Section (3.2). 

P and Q are measures of the errors in the constraints and the optimality conditions, 
respectively. For the exact optimal solution, one must have 


P = 0, 
Q = 0. 


( 9 ) 


For an approximation to the optimal solution, one must have 


P < 


(10a) 


Q < < 2, 


(106) 


where cj and f 2 are small pre-selected numbers. 

3.3.3 The Sequential Gradient Restoration Algorithm 

The Sequential Gradient Restoration algorithm involves a sequence of cycles, where 
each cycle has a gradient phase and may or may not have a restoration phase. 

The restoration phase is iterative and is started when the inequality (10a) is violated. 
In each restoration step, the norm of the variations of the control and the parameter 
is minimized, while the constraints are satisfied to first order; this has the effect of 
reducing P, i.e. more closely satisfying the feasibility conditions. The restoration 
phase is terminated when the inequality (10a) is satisfied. 

The gradient phase is started when the inequality (10a) is satisfied. It involves a 
single step, which reduces the functional 7, while the constraints are satisfied to first 
order; this has the effect of reducing Q, i.e. more closely satisfying the optimality 
conditions. 

The algorithm ends when both inequality (10a) and (10b) are satisfied. 
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Notation 

For steps in the gradient phase or the restoration phase, the following notation is 
used: x(r), u(r), tt denote the nominal functions; x(r), u(r), tt denote the varied 
functions, and Ax(r), Au(r), and A" denote perturbations of x(r), u(r) and 7r 
about the nominal values, i.e. 

i(r) = x(r) - Ax(t) , 

u(r) = u(r) + Au(r) , (11) 

7T = 7T -+■ ATT . 

(Note the removal of the arrow on top of the symbols for the vectors x, tJ and n in 
(11)) 

Concerning the functionals 7, J and P , the following terminology is used: 7, J and 
P denote the values associated with the nominal functions; 7, J, and P denote the 
values associated with the varied functions; and A7, A J, and A P denote the total 
variations of these functionals caused by the perturbations Ax(r), Au(r), and Att, 
i.e. . .. 

V- VC.- /= i t Ai , 

; - . : : .. . •; 

... .••:>*&« •• -#* 

; ■.. J=J + AJ, (12) 

P= P + AP. 

If the variations appearing in (11) are linear in the stepsize a , (where a > 0 ), they 
take the forms 


Ax(r)=Qw4(r) , 


Au(r)=aB(r) , 

(13) 

At r=oC . 

(14) 


3.3.4 Desired Properties in the SGR Process 


! 
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The functions A{~). B(r), and C are determined so as to decrease the functional I . 
and/or J, and/or P. (Alternatively, the condition (10b) on the value of Q can also 
be used.) Thus, the following properties are required. 

AI < 0, 

and or A J < 0. (15) 

and/or A P< 0. 


In turn, relations (15) can be enforced at every iteration providing the stepsize a is 
sufficiently small and the functions A(r). B(t). and C are chosen so that 



b I < 0, 

(16a) 

and/or 

b J < 0 , 

(166) 

and/or 

bP < 0, 

( 17 ) 


where 6(...) denotes the first variation. Inequalities (16a) and (16b) characterize 
the gradient phase, and inequality (17) characterizes the restoration phase. 

3.3.5 The First Variations 


The first variations of the functionals I , J, and P are 

bl = [ L t Ax + L u Au + L r An]dr + (<^ z Ar — <t> r A-n)\ 

Jo 


6J 


= [' \\ T + H x ]Axdr + f H u 

Jo Jo 

+ \(* z -\ T )Ax] l , 


Audr + j($,r)i 


[' H r dr) 

Jo 


At 


SP - 2 f (x - f) T (Ai - f x Ax - f u Av - f r A7i)dr 
Jo 

+ 2 [ S t (S z Ax + S u Au + SrAn)dT 

Jo 

+ 2[il’ T (^ x Ax + ^>rA»r)], , 


( 18 ) 


where (6x)o has been assumed zero in the expression for 6 J. 
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3.3.6 Gradient Phase 

Suppose that nominal functions x(r), u(r), and it satisfying (1) are available (a 
feasible solution). Let x(r), ti (r) , and it denote varied functions also satisfying (1). 
The varied functions are related to the nominal functions by (11). where Ar(r), 
Au(r), and A it denote the perturbations of x(r), u(r), and it about the nominal 
values. 


To first order, the perturbations Ax(r). Ai/(r). and Ait must satisfy linearized 
constraint relations 


Ax - f x Ax - /„ At/ - f K Ait = 0, 
S z Ax + S u Au + St Ait = 0 , 


(19) 


with the following linearized boundary conditions 

[Ax'o = 0 , 
jVxAx -+• xp t Ait)\ = 0. 


( 20 ) 


An inspection of equation (18) reveals that 6 J can be made negative through the 
following choice of variations of the control u and the parameter vector it. 


A u - -cxHl , 

Ait = -ai($„)i+ f Htdr} 7 . 

Jo 


( 21 ) 


where a denotes a scaling factor (gradient stepsize). The multipliers A(r), p[r), and 
H appearing in equation (21) must be consistent with the differential relation 

X T = -H x . (22 a) 


and the final conditions 


(A 7 '). = (*.)l • 


(226) 
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With the variations defined by (21). the first variation of the augmented functional 
(18) becomes 

SJ = -aQ. (23) 

where Q is the error in the optimality conditions (8), which, with (22). reduces to 

<3 = f' KtfjdT -r KU*lh ■> f Hi if, . (24) 

Jo Jo 


Since Q > 0, (23) shows that 6 J < 0. Hence, for sufficiently small a, the decrease 
in the augmented functional J is guaranteed. 


To simplify the problem, we can make use of the auxiliary variables A(r), B(t) 
and C introduced in (13). Using these variables, the linearized relations (19)-(23) 
become 


A — fzA — f u B — f*C — 0, 


(25) 


S X A t S U B ■+- S*C = 0 , 


with the following boundary conditions 

M)o = o, 
[tp z A + t/'irCji = 0. 


(26) 


The special variations are 


B= -Hi, 

C = -[(*,)! + [ Hxdr) 7 '. 
Jo 


(27) 


where once again the multipliers A(r), p(r), and n satisfy the following differential 
relation 

X T = ~H Z , (28a) 


and the terminal conditions 


(^)i = (*«)i. 


(286) 
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The descent property of the augmented functional (IS) is 

6 J = — a Q , 


where Q is the error in the optimality condition (24). which reduces to 


-i: 


B T Bdr - C T C 


(29) 


We note that the differential system (25)- (28) is linear and nonhomogeneous in the 
functions A(~), B(r), and C and the multipliers A(r), p(r), and p and can be solved 
without assigning a value to the gradient stepsize a. The selection of a is done a 
posteriori in such a way that the descent requirement (15-2) (second equation of 
equation (15)) is enforced. 

3.3.7 Solution Technique in the Gradient Phase 

The usual approach of backward integration of the A-equations and forward inte- 
gration of the .4-equations cannot be employed with this class of problems since the 
comptitftion of p(r) over the interval 0 < t < 1 requires the simultaneous solution 

: - : v y . y _ . 7i\y - 

of both the A and A equations. Because of this coupling, the total system must 
be integrated simultaneously either forward or backward. Forward integration is 
chosen here. 

Using the Method of Particular Solutions [39-40], the differential system consist- 
ing of equations (25), (26-1), (27-1) and (28) may be integrated forward n + p ■+■ 1 
times. Given 


A — f t A -l- f u B + f v C , 


( fl ) 


r\ — c a 


C D 


c n 

4 -' 7T , 




M)o = 0, 


(c) 


b = -hZ, 
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= -(£ W 

A 7 = -(£,+ A 7 /, - p T S,) . (e) 

Substituting equation (d) into (b) and solving for p. we get 

P = 'S u Slr'(S z A- S u Ll - S u fi\-S*C), (/) 

which is unique provided 

det S u Sj =£ 0. 

Therefore a necessary (but not sufficient) condition for the existence of a solution is 
that all the components of the path equality constraint vector S must involve some 
components of the control vector u. If this is not true, time derivative(s) of some of 
the constraints have to be taken to involve the control u in S. In principle, p may 
be eliminated by substituting (f) into (d), and (e) and B may then be eliminated 
by substituting (d) into (a). 

The initial conditions for the i th forward integration are 

(^)o = 0, 

A,(0) = \6 tl ,6 t2 6J t , ( g ) 

Ct — ^»(n + 2)’ •••’ ^t(n-t-p)] ’ (^0 

where i = l,2,...,[n + p + li and S tJ is the usual Kronecker delta function. Values of 
p and B can be computed using equations (f) and (g) respectively. The differential 
system: 

A = fxA ■+■ /uB + f n C , 

X = -L T t -fJ\-Sjp, 


can then be integrated forward. 
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The solution of the system is given by 

n-p- 1 

A(r)= £ M,(t). 

1=1 

n-p- 1 

B(r)= J2 MM')- 

1=1 

n-p- 1 

C = £ *« C »’ (30) 

1=1 

n-p-1 

*w= E m.m. 

1=1 

n-p- 1 

p( t ) = £ *!/>.(■)• 

1=1 

Here (...), is the value of (...) obtained from the i th iteration. The k t are unknown 
coefficients determined from the simultaneous solution of the following equations: 

n-p— 1 

= 1 , 

*=i 

n-p- 1 

) " ^i(V‘2-^t "" VjrC,)] = 0, 

1=1 

...... E M A ,). = (*Di • 

1=1 




where we have made use of the function 


$ = <P + /i 7 ^’ , 


defined earlier in equation (5). 

Equation (31) is equivalent to [l + g + n + p] scalar equations for the solution of 
the [n •+ p -i- lj unknowns k, and the q components of the multiplier //. Using these 
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coefficients from equation (31), the solution (30) satisfies the end constraints of (26- 
2), (27-2) and (28b). After the constants k, and the components of the multiplier p 
are known, the functions A(r). B(r), and C as well as the multipliers A(r) and p(r) 
are computed according to equation (30). In this way, the linear TPBVP is solved. 

3.3.8 Gradient Stepsize 

The gradient stepsize a is selected in such a way that the following inequality is 
satisfied. 

j (a) < j (0) , (32) 

subject to 

P(a)<f 3 , (33a) 

and 

t ; >0. (336) 

The last equation in (33b) expresses the need to have the unspecified end-time tf 
greater than zero. 

Here <3 is a small, preselected number (e.g. e 3 = 1 ). Satisfaction of (32) is guaran- 
teed by the descent property of the gradient phase. Satisfaction of (33a) is desirable 
in order to limit the constraint violation which is due to the use of the linearized 
constraint equations (25)-(26). Satisfaction of (33b) is automatic for problems with 
fixed terminal time, and is required in problems where the final time is free. 

3.3.9 Restoration Phase 

At the end of the gradient phase, the varied functions x(r), u(r), and t are known, 
and the varied constraint error P can be computed with equation (7). In this 
connection, two possibilities arise: either (i) Inequality (10a) is satisfied or (ii) 
Inequality ( 10 a) is violated. 
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Case (i) occurs if the equations (1) are linear or if they arc nonlinear and the 
gradient stepsize is sufficiently small. In this case, the gradient phase is repeated by 
employing as nominal functions x(r), t/(r), and n the varied functions x(r), fi(r). 
and ft of the previous gradient phase. 


Case (ii) occurs if the equations (l) are nonlinear and the stepsize is not small. 
In this case, prior to repeating the gradient phase, a restoration phase is inserted 
in such a way that (a) the constraint error is reduced to a level compatible with 
inequality (10a) and (b) the descent property of the augmented functional J is 
preserved. 

While the gradient phase involves a single step, the restoration is iterative and 
hence may involve several steps. This is due to the fact that the constraint equa- 
tions (1) are considered only in their linearized form in the restoration phase. For 
the first restoration phase, we employ as the nominal functions the varied functions 
x(r), «('), ft of the previous gradient phase. For any subsequent restoration it- 
eration. we employ as nominal functions the varied functions x(r), ti(r), ft of the 
pwvbus restoration iteration. 


Notv consider the generic restoration iteration. Let x(r), u(r), n denote nominal 

functions satisfying the initial condition x(0) = given, and let x(r), u(r), ft denote 

the varied functions also satisfy the initial condition. The varied functions are 

related to the nominal functions by (11). To first order, the perturbations Ax(r), 

Au(r), and Ait must satisfy the linearized constraint equations 
Ax - f t Ax - fuAxt - f n An + a(x - /) = 0 , 

(34) 


5* Ax + S u Au + S*An + aS = 0 , 


with the following linearized boundary conditions 

(Ax]o 


0 , 


[l£*Ai+ = 0, 


( 35 ) 
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where a denotes a scaling factor (restoration stepsize) in the range 0 < o < 1. Its 
function is to prevent the forcing terms in (34) -(35) from generating variations too 
large for the linearized assumptions to hold. 

When the variations defined by (34)- (35) are employed, the first variation of the 
constraint error (7) becomes (use equation (18-3)) 

6P = -2 aP, (36) 

where P is as defined in (7). Since P > 0. equation (36) shows that SP < 0. Hence, 
for a sufficiently small a, the decrease in the constraint error P is guaranteed. 

Since equations (34)- (35) admit an infinite number of solutions, the restoration 
iteration is not uniquely defined, unless some additional requirement is introduced. 
This added requirement is that the restoration be accomplished with the least square 
change in the controls and the parameters. Hence we seek the minimum of the 
following function 

K= f [Au T Au]dr + An 7 An, (37) 

Jo 

that also satisfies the linearized constraints (34)-(35). 

3.3.10 Special Variations 

The problem formulated in Section (3.3.9) is a linear-quadratic Bolza problem in the 
calculus of variations. The Au and At: that minimize (37) subject to the constraints 
(34)-(35) are: 

Au = -a(X T f u -+ p T S u ) T , 

r 1 ( 38 ) 

An = -a[ {X T U + P T S n )dr + {p T rpn)i} T ■ 

Jo 

The multipliers A(r), p(r), and p appearing in equation (38) must be consistent 
with the following differential relation 


X T = -(A T f x + p T S x ), 


(39a) 
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and the final conditions 

(A T )i = {p T ^x) 1 • 


(396) 


To simplify, we make use of the auxiliary variables .4(r), B(t) and C defined in 
(13). Using these variables, the linearized relations (34)-(35) become 


A-fzA-f u B-f w C~(z-f)= 0, 
S T A 4 S U B - 5*C 4 5=0. 


(40) 


with the following boundary conditions 

(A) 0 = 0, 

[WxA 4 1p*C 4 V 7 ] 1 = 0. 


(41) 


The special variations are 

B = — (A r / U -r p T S U ) T , 

C=-\f (\ T fn + p r S r )dT+ (Sxl’r)!? . 

Jo 


(42) 


the multipliers A(r), p(r), and p satisfy the following differential 

relation 


A r =-(A r / I = /5 I ), 


(43a) 


and the terminal conditions 

(A 7 ^ = {p T ipx)i - (^36) 

We note that the differential system (40)-(43) is linear and nonhomogeneous in the 
functions A(r), B(r), and C and the multipliers A(r), p(r), and p and can be solved 
wittmnt nceiunina a value to the gradient stepsize a. The selection of a is done 
a posteriori in such a way that the descent requirement (15-2) is enforced. This 
property is the same as in the gradient phase. 
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3.3.11 Solution Technique in the Restoration Phase 

Once again, because of the coupling due to the presence of p. the total system 
must be integrated simultaneously either forward or backward. Using the method 
of particular solutions, the differential system consisting of equations (40). (41-1). 
(42-1) and (43a) may be integrated forward \n - p + 1] times. Given 


A = ft A - f u B - / ff C - (i - /) , (a) 

0 = S Z A^ S U B ~ S T C + S , (6) 

(.4) 0 = 0 . (c) 

B — - (fy A -f- Sj p) . ( d ) 

A T = -(A T f I + p T S I ). (e) 

Substitution of equation (d) into (b) and solve for p. we get 

P = iS u Sjf 1 (S^-5 u /jA^S ff C^ S). (/) 


Note as in the gradient case that the above admits an unique solution providing: 

d et[S u Sy # 0 . 

Back substitution of (f) into (d) give 

£= -(/^A-Sjp), (?) 

where p is determined from equation (f). 

The initial conditions for the i*** forward integration are as given before (by equation 
(h) in the gradient phase). The differential system 


A = f x A + f u B + f„C - (i - /) , 
A = ~{ft A + Sjp) , 


(A) 

(0 


f 
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can then be integrated forward. The solution of the system is given by equation 
(30). The unknown coefficients k t , which differ from those in the gradient phase, 
are determined from the simultaneous solution of the following equations: 


n-!“p- 1 

L *. = i • 

f=l 

n+p~l 

+ ^-77^1)1 — ”(v)] , 

1=1 

1 

£ M^«) 1 = \v T {vi)]i • 

!=] 



sZt>,)ir + c,\ = -|((i 7 '(d>,)i! 7 \ 


(44) 


Equation (44) is equivalent to [l -f g ■¥ n p] scalar equations for the solution of 
the [r? -t- p — 1; unknown k, and the q components of the multiplier p. Using these 
coefficients from equation (31), the solution (30) satisfies the end constraints of (41- 
2), (42-2) and (43b). After the constants k, and the components of the multiplier p 
art known, the functions A(r). B(r), and C as well as the multipliers A(r) and p(r) 
art computed according ter equation (30). In this way, the linear TPBVP is solved. 

3.3.12 Restoration Stepsize 


The restoration stepsize is selected in such a way that the following inequality is 
satisfied. 

P{a) < P(0) , (45) 

subject to 


*/(<*)> 0. 


(46) 


In summary, the restoration stepsize must be chosen so that the above inequality 
constraints are satisfied. Should any violation occur, then a smaller value of a must 
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be employed and can be obtained, for example, with a bisection process, starting 
with q = 1. 

3.3.13 Summary of the SGR algorithm 

The SGR algorithm has cycles, composed of a gradient phase and (possibly) a 
restoration phase. The objective of each cycle is to decrease the functional J, while 
the differential and nondifferential constraints are satisfied to the predetermined 
accuracy required by (10). 

The gradient phase involves a single step, which reduces the functional J. It starts 
with nominal functions x (r), u(r), n which satisfy the constraints (1) within the 
preselected accuracy (10a). Using these nominal functions, the matrices f x , /„, /*; 
S z , 5„. S* etc. are evaluated at the integration intervals. The TPBVP (25)-(28) 
is then solved using the method of particular solutions. The functions A(t), B(t). 
C and the multipliers A(r), p(r), n, are used in turn to compute the corrections 
Ax(r), A u(r) and Atf after a suitable value of the stepsize a is found. The varied 
functions x(r), ti(r). n are given by equation (11). 

The restoration phase is iterative, hence involves one or more steps, and reduces the 
constraint errors to a level required by (10a). The nominal functions x(r), u(r), n 
are those from the previous gradient phase in the first restoration iteration, or those 
from the previous restoration iteration for any subsequent restoration iterations. In 
either case, they satisfy the given initial condition x(0). 

Using these nominal functions, the matrices f z . /„. /*; S z , S tt , S, r , as well as the 
vectors S and (i — /’, are evaluated at the integration intervals. The TPBVP (40)- 
(43) is then solved using the method of particular solutions, obtaining the functions 
w4(r), B(r), C and the multipliers A(r), p(r), and /*, which are used in turn to 
compute the corrections Ax(r), Au(r), An after a suitable stepsize a has been 
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found. The varied functions x(r). u(r), ir can then be computed using equation 

( 11 ). 

The satisfaction of the condition (10b) signals the termination of the process and 
the converged solution is given by the varied function x(r), u(r), n at the end of 
the restoration phase. 

A flow chart of the SGR algorithm is given in Figure (3.3.3). 

3.3.14 Example Problem 


The example given here is a modification of the autorotation problem that was 
formulated in Section (2.9). Here we add an upper bound on the vertical sink rate. 
Experience gained from the study on the optimal landing of a helicopter after engine 
failure indicates that the optimal control program usually results in an unacceptably 
high maximum sink rate. This upper bound on sink rate appears as a state variable 
inequality' constraint: 



v deicent — {^descent) 


maximum ? 


(47) 


which can be transformed to a path equality constraint by the use of a Valentine- 
type slack variable. Since i\ denotes the normalized vertical sink rate, we have 


*i ~ Ilm + *6 = 


(48) 


where xg is an auxiliary state variable created to convert (47) to (48). Note that 
the initial condition of X6 is given by 


Te ( 0 \ = 


• */h 

- v r 


-f'fnll 


(4Rn) 

\ ~ ~ ' / 


where either the “-j-” or " sign can be used. 
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(3) Gradient restoration cycle 

(4) Stability cycle 

Figure 3.3.3 Flow Chart of the Sequential Gradient 
Restoration Technique 
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The first time derivative of (48) is 

Xi' + 2i 6 t/4 = 0, (49) 

where ( )' denotes time differentiation with respect to the normalized time r. and 

u 4 is an auxiliary control variable defined as 

= Xg. (50) 

Since 

x'j = go - m 0 (t/ix^ + f x \\' x \ 2 + * 2 2 ) , (51) 

we can substitute (51) into (49) and convert it into the following path equality 
constraint that involves the two controls uj and u 4 : 

9o - TO 0 (uix 2 + f x i V 1 ! 2 + X 2 2 ) + 2x 6 ti4 = 0 . (52) 


The new form of the autorotation optimization problem with vertical sink rate 
bemad is as given below. Let 


X = (xj X2 X3 X 4 X5 Xo) T , ( 53 ) 

U = (ui u 2 u 3 u«) , (54) 


7T = (r ; ). 


(55) 


denote the augmented state, control, and unknown parameter vectors. 

The problem is to find U{£) and f to minimize: 

I = \{ x \* + w z x 2f 2 ) = <t>(X f ) , (56) 


subject to : 
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(1) equations of motion ( A v = f ): 



II 

o 

1 

3 

o 

p» 

(]X 3 2 - /x 1\ Xi 2 - x 2 2 )), 


= Tjmo{u 2 X 2 2 ■ 

- /x 2 \ X ) 2 -r X 2 2 ), 


= -Tfioxz 2 {co 

- A \/v 1 2 - t/ 2 2 ), 

* 4 V 

= 0 . 1 7 f X ] , 


Ih V 

= 0.177x2. 


*c, V 

= TfU 4 . 



(2) the initial condition: 


Xq = (0, u 0 , 1,0,0, v /i lm ) 7 '. 


(58) 


(3) path equality constraints (S(A\ £/,?:) = 0) are: 




2 (uj 2 -f t/2 2 )v 2 


V-V 3 2 = 0 ., 


Ct. 


go ~ m 0 (u] x 2 + fx\ \'xi 2 4 - x 2 2 ) - 2x 6 U4 = 0. 


(4) terminal constraints (v'(A”/,7r) = 0) are: 


x Af - hf = 0, 

xsf ~ df = 0 . 


(59) 


(60) 


Simplification 

Since the computational effort involved in obtaining the optimal control solution 
using the SGR algorithm is approximately proportional to the square of the dimen- 
sion of the state vector, the optimal landing problem in its “reduced” form is much 
preferred. 
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Since equation (48) is linear in the state x\. we can find x\ in terms of the other 
state variables and eliminate it from equations (53)-(60). We first note from (48) 
that 

T\=x\m-A- ( 61 ) 

Substituting (61) into equations (53)-(60) gives the following reduced problem. Let: 

X = (l 2 *3 *4 *5 a*) 7 , ( 62 ) 

U = (l/] V2 ti3 U4) r . (62) 

7T = (Tf). (64) 

denote the reduced state, control, and unknown parameter vectors of the problem. 


The problem is to find U{£) and ir to minimize: 


/ = 5l(x lm - x 6 5) ! + W'.x,/] = t(X ,) , 


2 ] _ 


(65) 


subject to : 

(l) equations of motion ( A^ = / ): 


X2 V 

= r/mo(«2X3 2 

~ f x 2\/(Xlm - x\ 

X3 V 

= — r/*oX3 2 (co 

+ X\/ Ui 2 -r U 2 2 ), 

X4 V 

= 0.1r / (i lm - 

*}). 

X5 V 

= 0.1r^i2, 


xe V 

= TfU A . 



( 66 ) 


The equation of motion for x \ , which is not needed in the reduced formulation has 
been removed. Note also that the inflow ratio A used in (66) is given by a modified 
version of equation (43) (cf. Chapter 2) given below: 


A = O-Oli 1 ’" 8 }*‘ m J ; )Ul ] + *o//(*i 2 + «2*)i . 

X3V U 1 2 + u 2 2 


( 67 ) 



Results obtained for this problem are given in Chapter 4. 

§3.4 Neighboring Extremal with Path Constraints 

The necessary conditions for optimal control problems with path equality con- 
straints are given by equations (4)-(6) of the last section. The SGR algorithm 
described in that section is one way that such a problem can be solved. However, 
this solution is open-loop and even slight changes in the initial or final conditions 
require a new computation of the entire control program. 

A linear neighboring optimal feedback control scheme was developed by Breakwell, 
Speyer, and Bryson for optimal control problems in which the state variables are 
subjected to initial and terminal constraints [50’. Their approach minimizes the 
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second variation of the performance index. Kelley developed a similar feedback 
control scheme using the second variation in an accessory minimization formation 
(51’. Bryson and Ho [29 utilized a Riccati transformation with a backward sweep 
method to obtain the feedback gains needed to approximate neighboring optimal 
solution. 

In this sub-section, we extend these earlier works to problems which include path 
equality constraints, and specifically to problems with both path equality con- 
straints 5 and an unknown parameter vector if. The necessary conditions for this 
class of problems are given by equations (4)-(6) of Section (3.3). These necessary 
conditions are repeated here for easy reference: 

X = /(£,U.7f,7-) . (1) 

S(x,tT, ir,r) - 0, 

■ T 

A =-H z , 

H u = 0, 

x(0) = given , (2) 

(A 7 ), = (*«)i. 

^(x;,7f) = 0, 

and ($*)i *r f H n dr = 0. 

Jo 

Now consider small perturbations from the nominal optimal path produced by small 
perturbations in both 6x(0) and 6xp. These perturbations give rise to perturbations 
6x(r), 6\(t), 6u(r), 6p[r), Sfi and Sir. The relationships among these perturbed 
quantities are given by 
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(l) Perturbed differential relations: 



(2) Perturbed path equality relations: 



o = ($,ri)] (6x)i - (tV r )i(£/7) - : ($»T7r)i f H nr dr]67r 

J 0 

-(/ H T1 6xdr)-([ H TU Sv dr) - ( f fj6Xdr)~(f Sjbpdr) . 
Jo Jo Jo Jo 


(3) End conditions 

6f(0) = specified, 

(6f)i = (t,)i(6f)i - (vxjjtf = specified. (6) 

(«A), = (*„)i(«i)i - + (v-JhSp. 
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From equation (4), we can solve bu and 6 p in terms of 6 x, 6 A. and Sir as follows: 



A\2 

A 22 


A 13 
-*23 



(7) 


Equation (7) can be substituted into equations (3) and (5) and the result is the 
following linear, time-varying differential equation 


6 x 

bX 


C\ 1 C \ 2 C \3 

Ci\ C22 Ciz 



(8) 


with the integral path constraint 


/: 


5jo + SnSn / (S]2^x-h S\ 3 6 A)<fr = 0 . 


(9) 


Here A tJ , C tJ (j = 1.2 and j = 1,2.3) and 5io-5i3 (not to be confused with the path 
equality constraints 5) are time- varying matrices evaluated on the nominal optimal 
path. Detailed expressions of these matrices, as well as those of H uu , H ur , H zz , 
HuxiHuz and H ru are given in Appendix (E). Using the results given in Appendix 

t'.j. . • • • '• , v 

■■■■■•'■ T 

(E), we can easily establish that -C\\ = C22 and that both C\2 and C21 are 

symmetrical matrices. 

Note that the system of equations (3)-(6) admits an unique solution providing: (l) 
H uu is non-singular; and (2) S 4 (r x r) (= S u is also non-singular. 

The first requirement is the usual convexity condition (or strengthened Legendre- 
Clebsch condition). The condition is easily understood from the fact that u(r) is 
determined by minimizing the Hamiltonian function ( H ) with respect to u while 
holding x, A and if fixed. For a smooth H- function with no constraint on u(r), it 
requires that 

H u = 0, and H uu > 0. 
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The second requirement ( i.e . S 4 > 0) can thus be interpreted as a modified form 
of the convexity condition for problems with path equality constraints S. Note 
that if any of the path equality constraints does not involve some components of 
the control vector t7, the determinant of S U H~J Sj i is zero. A necessary (but not 
sufficient) condition for the non-singularity of the S 4 matrix is hence the involvment 
of the control u in 5. 

Simplifications 

If <?. t\ S and f„ are not functions of ft. some simplifications are possible. These 
include: 

(1) Simplified expressions for some of the 5-matrices (cf. Appendix (E)): 

■^2 = S U H UU H U jr , 

Sj = 5jo = 0 , 

Sn= f ( H UT S 6 )dr , 

J 0 

S 12 = Hm ■+■ Hn U Sf, , 
and 5 j 3 = H nu S: - /J . 

(2) Simplified end-conditions: 

6f(0) = given , 

(ftp) 1 = ( Ci ) ] ( 6 x ) 1 = given , 

(<A),= (*„)i(M)i + (<£),(*;<)• 

(3) Simplified expressions for H^y 

Since H(x, u, ft, A, p, r) - L(x, u, 7r, r) -f A T f(x, o, ft, r) - 1 - p T S(x, u, r) , 
therefore H u = L u -t- A T f u -|- p T S u , 

H x = L x -i- A t fx ■+ (? S z , 

H r = Ijr + A 7 /, . 
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General expressions of H ( . ) are given by 

H.,= |Vj) = + 

B.. = ^(Hl ) = ^(Ll + ffX- Sjo) , 

*..= = ^(iJ + ArA). 

»«. = = ^Ll + S*p ) , 

if ’ ,= 

and H„„ = ^(ii T + /, T A). 

With these simplifications, the linear TPBVP (given by equations (6)-(9)) is given 
in Section (3.4.1). 

3.4.1 Problem Formulation 



with the following integral path constraints 

SnSir -r / (S 12 6 x + Si 3 6A)dr = 0, 

VO 

and subjected to the following end conditions 

6£(0) = given , 

(Stf) i = (0*)i(6x)i = 5J'ven, 

(«A)i = (*«)i(*f)i + (lk J V*. 

3.4.2 Solution by the Backward Sweep Method 

6\ and 6xp can be expressed in terms of 6x , 6jl and 6n in the following form (cf. 
Chapter (6) of reference [29]) 

6 X= T6x + Roii + Foil, 

6 ip = R T 6 x + Q 6ji + H 6i F . 


( 11 ) 
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Here T(n * n), R(n y q), P(r> > p). Q(q v 9) and #(9 • p) are matrices that satisfy 
the following differential equations: 

f = -TC,i - Cj.T-TCnT^Cn . 

i?= -(cS-rTCnJi?, 

P= C 23 - Cf, P - T(C n P * C, 3 ) , ( 12 ) 

C? = -f? r C 12 i?, 
and H = - R T {C n P - C, 3 ). 

The first three relations in equation (12) can be obtained as follows. We first take 
the time differentiation of 6 A in equation (11) and substitute into the resultant 
equation expressions of 6 i and 6 \ from equation (10) (noting that both Sjl and 6 t: 
are constant). The final equation has the form 

[ ] 6 x — ! ] 6 jl — ! i^ 7 r = 0. 


If the above equation is to be identically true for all Si. Sp and Sir, all the ’ 1 terms 
must be zero. Thus, we obtain the first three relations of equation ( 12 ). 

The last two relations in equation ( 12 ) are obtained in a similar way. Here we note 
that Stp is a constant when we take the time derivative of the second equation in 
(11). Substitution into the resultant equation expression for Si from equation (10) 
and the expression for SX from equation (11) give the differential relations for Q 
and H. 


The terminal conditions of the T, R, P, Q and H matrices are obtained in this 
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manner and are given by the following relations: 

r(i)= 

(vDi, 

P(l) = (*«)i = 0, (13) 

Q( 1 ) = 0 , 

and H( 1 ) = (v»)i = 0 . 

Finally, to obtain the neighboring optimal feedback law, we can substitute equation 
(11) into equation (7), yielding the following relation: 

6 U = — K\6x — A-iSrp — A 367 F , (14) 

The time- varying neighboring optimal feedback gains, A, ( t = 1, 2 and 3) are given 
by the following relations: 

Ai(mxn)= -\An + An{T - RQ~ l R T ), , 

A2 (m x 9) = -Ai2RQ ~ } , ( 13 ) 

V, .7- *nd A 3 (m x p) = -\Au + An{P - RQ~ X H)\ . 

Equation (14) is a continuous linear feedback law that will produce the desired small 
changes in the terminal conditions 6tl> starting from the revised initial conditions 
of x( 0 ) -f 6 f( 0 ). These end-conditions are met in such a way that the performance 
index 1 (and hence J) is minimized. This method constitutes the Neighboring 
Optimal Feedback Law. 

In addition to our earlier mention of the non-singularities of the H U u and S U H "JSj - 
matrices, we now note from equation (15) other necessary conditions for a local 
minimum in J. These additional necessary conditions are that both 


\T-RQ~'R t \ 
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and P- BQ-'H 

be finite for 0 < r < 1. This is similar to the classical Jacobi condition, or the 
condition that no conjugate points exist on the path. If either of these conditions 
is violated, for example, ; T - RQ~ l R T ; approaches infinity at time t' (where 0 < 
t' < 1), it is necessary that certain combinations of the components of b x at time r' 
be zero. The set of all possible perturbations is now restricted to a dimension less 
than n (the dimension off). In fact, the path is not a minimum path if it continues 
beyond r'. 

Note from equation (13) that the matrix Q is identically zero at the terminal time. 
All the gain matrices given by equation (13) approach infinity as the end-time 
nears. This mathematical singularity can be bypassed if we adopt the strategy of 
“switching-off” the feedback law’ as the end-time is approached. 

The value of df is needed in the neighboring optimal feedback law. This value must 
be determined by the simultaneous solution of both the linear TPBVP and path 
equality constraints. Appendix (F) gives an iterative procedure that can be used to 
find bit by taking advantage of the fact that bit enters both the differential equation 
(10) and its initial condition linearly. A schematic diagram of the neighboring 
feedback control is given in Figure (3.4.4). 

3.4.3 Example Problems 

Two example problems are given here to demonstrate the concept. The first is for 
a dynamic optimization problem constrained by a path equality constraint 5, but 
without the variable parameter vector tt. The second example includes both 5 and 
H. The computer code that is used in the following computation is called Program 
“SECOND”. Furhter information about the computer program can be found in 
Appendix (F). 



Figure 3.4.4 Neighboring Optimal Feedback Law 
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Example Problem (1) 


Problem Statement 

min 1 = —i jf 
e 1 

x = V cos 6 , 
y = V sin 6 . 

The initial conditions of the problem are 

*( 0 ) = 0 , 
y( o) = o, 


and the terminal constraint 

x(t f ) = x f . 

The end-time of the problem tj is fixed. 

This is a simple two-dimensional problem of transferring between the initial and 
final states of a vehicle with a constant velocity V. The objective is to control the 
velocity orientation angle 6 in such a way as to maximize y(t/). The problem is 
simple enough that analytical solutions for both the optimal control problem and 
its neighboring optimal feedback law can be determined as shown below. 


(A) Optimal Control Solution 
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Using standard optimal control techniques (cf. Chapter (2) of [29 ). the optimal 
solution is given by the following relations: 

sin*' = (l-liL)’)*, 

‘ </ 

cosr = (^) 
r’(') = 

u 

»*(0 = v(i-($-)*)h. 

Therefore the optimal orientation of the velocity vector 6* is fixed and is given by 
the above expression. The optimal trajectory' in space is a straight line and the 
problem is well posed only when 


I xj !< vtf. 


(B) Neighboring Optimal Feedback Control 

The neighboring optimal feedback control law can also be determined in the usual 
way (cf. Chapter (6) of (29j). The optimal feedback law is given by the following 
equation: 

66' = -Ai 6x - AtSip . 

Here Aj(l x 2) and A 2 (1 * 1) ere the feedback gain matrices associated with the 
state and the terminal constraints respectively. These gain matrices are given by 
the following relations: 

A i = ( v(i/- °) ♦ 

A 2 = ( - ) • 
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Substitutions of these gain matrices into the feedback law give 


£ 0 . _ _/ ^ \ \ 

[ Y{lj - t)sin6* * (// - <) sin 6‘ 

These “exact" analytical results are used as references for comparison with numer- 
ical results. 

The above problem can be re-formulated into one with path equality constraint as 
follow 


Alternative Problem Formulation 


min I = — xi{t f) 

€ 

i l = 2 ui , 

X2 — 2u2 . 

The initial conditions of the reformulated problem are the same as those given 
before 

*i(0) = 0, 

*2(0) = o, 

The terminal constraint is 

X(tf) = 1. 


The new problem is the same as the original one with V = 2, tj = 1 and xj = 1. 
Instead of using the orientation angle 6 as the control, we choose to use the sine 
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and cosine of the same angle as the controls of the problem. That is, 

U] = COS0, 

t/ 2 = sin# • 

Because of these control definitions, the problem is now constrained by a path 
equality constraint 

Uj •+■ v% = 1 . 


The neighboring optimal feedback law of the problem in its alternative form is solved 
using the numerical technique described earlier. The procedure is started with the 
solution of the nominal optimization problem and followed by the computations of 
the neighboring feedback gain matrices Aj and A 2 . The non-zero components of 
these matrices are given in Figure (3.4.5). 


These numerically obtained gain matrices should now be compared with those de- 
termined analytically. With if = l.V = 2, and tj = 1, we have from the analytical 


results 


sin 

or cos 6‘ 


\'b 

<!>• 


and S6‘ = -( 


6 x 


) + (- 


bit 


>/3(i-t) V5(i-<) 


) 


But since 

Uj = cos 0, 
U 2 = sin0. 

therefore we have 
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Substitution into the last equation expression for 66 * gives us the neighboring op- 
timal feedback law, (6uj . 6 u 2 ) T 



These results agree with those obtained numerically (see also Figure (3.4.5)). Note 
that the magnitudes of these gains approach infinity at the terminal time. 

The usefulness of the neighboring optimal feedback law can be confirmed in sev- 
eral ways. For example, the initial and terminal conditions of the nominal control 
problem may be “perturbed” simultaneously as follow 

xj(0) = 0.0— +0.1 ; 6xi(0) = +0.1, 

X2 (0) = 0.0— —0.1 ; 6x2(0) = — 0 . 1 , 

Xj (l) = 1.0— +0.9; 6xj(l) = -0.1. 

Th* optimal iu>ltttioo to this perturbed problem can be obtained numerically using 

SECOND. The solution may also be determined analytically using the revised end- 
conditions. 

Figure (3.4.6) shows the results obtained numerically for both the nominal and the 
perturbed problems. The former is obtained using the sequential gradient restora- 
tion (SGR) technique while the later by SECOND. It can be seen from these figures 
that the neighboring optimal feedback law has predicted substantial changes in both 
ui and u 2 , resulting in large changes in the time history of X](t) and X2(0- 

Both the analytical and numerically computed results for the perturbed problem 
are given in Figure (3.4.7). The figure shown indicates excellent agreement. 


Example Problem (2) 



Nominal Path 

: Predictions of Neighboring Optimal Fmmdbaek Laws 


Figure 3.4.6 Nominal and Perturbed Results in Example (1) 








Direct Optimization Technique with revised and-conditions 
Pradictions of Neighboring Optiaal Faadback Laws 


Figure 3.4.7 Comparison of Exact Results with those from 

w 11 tr •_ T* 1 _ / « \ 

reeaoacK Law iu 
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The second example is one with the terminal time unspecified. The problem state- 
ment is 

min 1=1 \ u \ ^ x \ dl 

u Jo 

i\ = ui , 

X 2 = («2 — Xj ) . 

The initial conditions of the problem are 

X) (0) = 0, x 2 (0) = 0. 

The path constraint is given by the following equation 

u\ + Xj = 1 . 


and the terminal constraints are 

x,(f/)=0.5, * 2 (</)= 0.5. 

Note that the problem is linear except for the path constraint. The problem is 
designed in such a way that both H uu and do not vanish for 0 < t < tj. 

x(2 x 1). u(2 ✓ 1), t.’(2 >- 1) and S(1 x 1) are respectively the state, control, terminal, 
and path equality constraint vectors of the problem. 

The above problem can be transformed into one with fixed end-time with the fol- 
lowing change of independent variable 

t 

T ~ h 

The independent variable of the problem is now r which varies from 0 to 1. 

The transformed problem is 

min/ = / tt\u\ + x\\dr 

« Jo 
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X\ = tfV i , 

*2 = tf(u 2 - Xj) . 

The initial conditions of the problem are 

xi(0) = 0, x 2 (0) = 0. 

Here ( )' denotes differentiation with respect to the new independent variable r. 
The path constraint is given by the following equation 

«2 + *i = 1, 


and the terminal constraints are 

xj(l) = 0.5, x 2 (1)=0.5. 

We have introduced in the above process an unknown parameter t j into the problem. 
The only component of the vector if is hence the unspecified end-time of the problem 

</• 

The nominal optimization problem is first solved and the neighboring optimal feed- 
back gain matrices (Aj , A 2 , A3) in the following feedback law are computed. 

Su(r) = -Ai(r)6x — A(r)6i/> - K.$(t)6t . 

Here Aj, A 2 and A3 are (2 x 2), (2 x 2) and (2 x 1) matrices respectively. The 
non-zero components of these matrices are given in Figure ( 3 . 4 . 8 ). Note that (once 
again) these gain components approach infinity near the end-time of the problem. 


The nominal optimal control problem is being perturbed in the following ways 
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Figure 3.4.8 Feedback Gains in Example (2) 
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(1) Perturbed initial conditions 

6x[0) = (0.05,0.05) r 

(2) Perturbed terminal conditions 

<5v(l) = (0.05, -0.05) 7 

(3) Perturbed initial and terminal conditions 

6x(0) = (-0.03, 0.03) 7 , 6v~(l) = (0,03, -0.03) 7 . 

Since analytical solution of the problem cannot be easily obtained, these problems 
with revised end-conditions must be solved by using the SGR technique. The results 
obtained in this way are considered to be “exact” . The neighboring optimal feedback 
law is then used to solve the same problems. In this way, the accuracy of the results 
obtained by using the feedback law can be gauged by a comparison with the “exact” 
results. 

(1) Per<urb«d l»ftfa] conditions 

The time history of both x(r) and u(r) in the nominal control problem are given in 
Figure (3.4.9). Those for the perturbed problem with revised initial conditions, as 
predicted by the feedback law are given in the same plots. The purpose of showing 
these figures is to give some indication of the percentage changes in the state/control 
vectors that have occurred due to the revision of the initial conditions. 

By using the SGR technique, the end-times of the nominal and perturbed problems 
are found to be 0.728 and 0.681 respectively. Therefore the change in the terminal 
time t f is 
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ORIGINAL PASuJS 

OF POOR QUALITY 






Predictions of Neighboring Optimal Feedback Laws 


Figure 3.4.0 Nominal and Perturbed (initial condition) Results 
in Example (2) 
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Using the iterative procedure described in Appendix (F), the value of btj is found 
to be —0.045 (with t = 10 -6 , see Appendix (F) for definition of f). The computed 
bt f is in good agreeement with the “exact” value. Error between these values of 
btj is mainly due to the computational error incurred in the backward integration 
of equation (12) and also due to the early termination of the iterative procedure at 
£ = 10 ~ 2 * * * 6 . 

The solutions to the perturbed problem obtained by using the neighboring optimal 
feedback law and those obtained directly by using the SGR teennique are plotted 
together in Figure (3.4.10). These plots shown excellent agreement between results 
obtained from the two approaches. 


(2) Perturbed terminal conditions 

The perturbations of the terminal conditions are by = (0.05 , —0.05)^. The neigh- 
boring optimal feedback law is once again used to compute the revised optimal 
control program. By using the SGR technique, the end-times of the nominal and 
perturbed problems are found to be 0.728 and 0.700 respectively. Therefore the 
value of btf is given by 

bt f = 0.700-0.728 = -0.028. 

Once again the iterative procedure of Appendix (F) is used to estimate the change 
in the end-time bt f and the value computed is -0.0278 (with £ = 10 -6 ). This value 
agrees well with the “true’’ value of —0.028. 

Since Ai, A2 and A3 approach infinity near the end-time, we adopt the strategy 
that whenever bu(r ) is larger than a pre- determined amount, the neighboring op- 


> 



Direct Optimization Tachniqua with revised and 'conditions 
Predictions of Neighboring Optimal Feedback Laws 

fit f (Direct) « *0.047 
fit t (Predieted) • *0.045 


Figure 3.4.10 Comparison of Exact Results with those from 
Feedback Law in Example (2) 
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timal feedback law should be turned “off". Figure (3.4.11) shows the time history 
of x(r) and u(7) of the nominal control problem and those of the perturbed prob- 
lem as predicted by the neighboring optimal feedback law. The usefulness of the 
neighboring optimal feedback law can once again be confirmed by comparison of 
the results obtained by these different approaches. Such a comparison is given in 
Figure (3.4.12). 

Other than the time history of ui(r), the results given in Figure (3.4.12) show 
good overall agreement with the “exact” results. The prediction for 6ui(r) is very 
good initially but deteriorates as the end-time is near. This error near the end- 
time is quite common among neighboring extremal problems with revised terminal 
conditions. In the present case, the difficulty is compounded by the the fact that 
St j has been “inaccurately” determined (Stj (estimated) = —0.0278 as compare with 
-0.0280). This small error (-0.0002) in the determination of Stj is being amplified 
by the magnitude of As( t/j component ) which becomes very large as the end-time 
is approached. Therefore the neighboring optimal feedback law over-estimates the 
value of u\ by a substantial amount near the end-time. Since «2 component ) is 
identically zero, the catenation of «2 is not being affected by this error in Stj. 

(3) Problem with perturbed initial and terminal conditions 

Here the initial conditions are perturbed by Sx{ 0) = (-0.03,0.03) r while the ter- 
minal conditions are perturbed by Stp = (0.03, -0.03) r . The neighboring optimal 
feedback law is used to obtain the optimal program for the perturbed problem. 
Using the SGR tecnhique, the end-times of the nominal and perturbed problems 
are found to be 0.728 and 0.650 respectively. Therefore the value of 6tf is given by 

Stj = 0.650 - 0.728 - -0.078 . 

Using the iterative procedure of Appendix (F), the change in the end-times of 
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Nominal Path 

Prediction* of Neighboring Optimal Feedback Lava 


Figure 3.4.11 Nominal and Perturbed (terminal condition) Results 
in Example (2) 
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Direct Optimization Technique with revised end-conditions 
Predictions of Neighboring Optimal Feedback Laws 

, (Direct) » -0 0280 
6t , (Predicted) « -0 0278 


Figure 3.4.12 Comparison of Exact Results with those from 
Feedback Law in Example (2) 
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the nominal and perturbed problems is computed to be —0.077 (with < = 10 -6 ). 
The estimation of 6i j is once again good. Figure (3.4.13) show's the time history 
for both the nominal control problem and the solution to the perturbed problem 
obtained using the neighboring optimal feedback law. Figure (3.4.14) shows the 
results obtained for the perturbed problem using either the neighboring optimal 
feedback law or the SGR technique. The overall agreement between these results 
is good. The cause of the discrepancy found in the time history of uj(r) was given 
before. 

3.4.4 Conclusions 

Neighboring extremal paths for dynamic optimization problems with path equality 
constraints and unknown parameter vector were considered in this section. With 
certain simplifications, the problem can be reduced to one of solving a linear TPBVP 
with integral path equality constraints. The solution technique employed is an 
extension of the backward sweep method given in Bryson and Ho (cf. reference 
[29). 

Two example problems are solved using a numerical algorithm called SECOND. Ex- 
cellent agreement with the “exact'’ results can be obtained for problems with small 
changes in the initial conditions. For cases with changes in the terminal conditions, 
the agreement deteriorates when the end-time is approached. This problem near 
the end-time is compounded by the multiplication of the small error involved in the 
determination of 6W and the large feedback gain (A,) near the end-time. One prob- 
able way of overcoming this problem is to switch from problems wdth hard terminal 
constraints into problems with only soft constraints as the end-time is near. This 
idea has not been investigated in depth and is recommended for future research. 
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Noainal Path 

Predictions of Neighboring Optimal Feedback Laws 


Figure 3.4.13 Nominal and Perturbed (initial and terminal conditions) 
Results in Example (2) 








Direct Optimization Technique with revised end-conditlone 
Predictions of Neighboring Optimal Feedback Laws 


6t f (Direct) « *0.077 
fit r (Predicted) ■ *0.076 


Figure 3.4.14 Comparison of Exact Results with those from 
Feedback Law in Example (2) 






Chapter 4 

Optimal Solutions and 
their Interpretations 


In this chapter, optimal solutions are obtained for the the helicopter autorotation 

problem without and with the descent velocity bound (cf. Sections (2.9) and (3.3)) 
using the algorithm developed in Chapter 3. 

We begin with a description of an OH-58A helicopter modified with a High Energy 
Rotor System (HERS). An autorotation flight test program [12 was conducted 
by the Bell Helicopter Company using this helicopter in 1976 '77. The data from 
this ;|e^gram ate suitable for comparison with analytical results of this research. 

• »* •• . .•‘V- •> v.Vi • - " « - 

Sections {4.2} and (4.3) describe the energy management considerations during au- 
torotation, and the piloting techniques that were employed during the autorotation 
maneuver. These techniques are compared w’ith those associated with the optimal 
control schemes. 

Optimal solutions obtained for autorotation landings from hover, from forward 
flight, and from forward flight with a descent velocity bound are given in Sec- 
tions ( 4 . 4 ), ( 4 . 5 ), and ( 4 . 6 ) respectively. These results are presented using plots 
of the time histories of the thrust coefficient, collective pitch, and states of the 
helicopter. In each case, the optimal results are interpreted physically and the ad- 
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vantages and disadvantages of the optimal control scheme are discussed. We end 
with some general comments on the effectiveness of the algorithms used in solving 
dynamic optimization problems with path inequality constraints. 

§4.1 Description of Test Vehicle 

The helicopter model used in our study of autorotative landings is a point-mass 
model of a modified OH-58A, shown in Figure (4.1.1). This helicopter was equipped 
with an experimental High Energy Rotor System (HERS), and was constructed by 
the Bell Helicopter Textron Company (BHT) as a concept demonstrator. The 
inertia of the rotor system was increased by the addition of tip weights in the spar 
cavity of the blades (hence the name “High Energy Rotor”). Flight evaluations [12 
indicated that the high energy rotor could eliminate the usual height-velocity (H- 
V) restriction, and that the rotor kinetic energy could be used to provide transient 
power for better maneuverability. This model is used in our study because flight 
test data of the helicopter in autorotation are available in references [11], [12] and 
[57]. 

The HERS consists of a two-bladed rotor with a diameter of 35.3 feet which operates 
at 354 rpm. The rotor blades have 16 inch chords and were designed so that the 
rotational inertia of the rotor system could be varied from the standard OH-58A 
inertia to over twice this value. Four external steel doublers were added to carry 
the increased centrifugal force. Lock numbers, of 5.43, 3.19, and 2.61, which 
correspond to blade inertias If, of 323, 550, and 672 slug- ft 2 respectively, were used 
in flight tests. The Lock number, of the rotor is defined as 

pacR 4 

which is a measure of the ratio of aerodynamic and inertia forces on the blade. 
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Here 

p = density of air (slug/ft J ) , 
a — lift curve slope of rotor airfoil (rad 1 ) , 

c = chord of blades (ft) , (2) 

R = rotor radius (ft) . 

/(, = rotational inertia of one blade and hub (slug ft 2 ) . 

Another quantity, the rotor solidity a, is of importance to the autorotational per- 
formance of helicopter. The rotor solidity is defined as the ratio of the total blade 
area (NcR for constant chord blades) to the total disk area (7r/? 2 ). Therefore, the 
rotor solidity is given by the following expression: 

Nc 


c = 


n R ' 


( 3 ) 


where ,V is the total number of blades per rotor. The solidity of the HERS rotor is 
0.048. 


A detailed description of the helicopter is given in reference [56 r Table (4.4.1) gives 
the values of the model parameters used in the point-mass model of this vehicle 
that was developed in Chapter 2. 

4.1.1 Height-Velocity Restriction Curves 

Most helicopters have a region of operation from which a safe autorotational land- 
ing cannot be executed. This area of limited autorotational capability exists for 
both single and multiple engine helicopters and is described by the altitude above 
the ground and the airspeed. It is commonly illustrated with a height-velocity re- 
striction diagram (sometimes refered to as a deadman’s curve). The height velocity 
restriction curves for the standard OH-58A helicopter, as determined by flight tests 
and reported in reference [57], are given in Figure (4.1.2). The two restriction re- 
gions indicated in Figure (4.1.2) are typical for conventional helicopters. The low 



S/NO. 


SYSTEM PARAMETER 


VALUE USED 


1 

V equivalent flat plate area (ft 2 ) 

24.0 

2 

g, acceleration due to gravity (ft/eec 2 ) 

32.17 

3 

P. air density (slug/ft 3 ) 

0.002378 

4 

a. aass of helicopter (slug) 

(lb.vt.) 

93.16 

(3000.0) 

S 

a. rotor solidity 

0.048 

6 

K , induced velocity correction factor 

1.13 

7 

*• £< c < > 

5.73 

8 

R, radius of main rotor (ft) 

17.63 

;rr- 

* '.t^v 

• ' *>. . 

Ji speed (rpa) 

354.0 

io 

6.. profile drag coefficient 
(NAGA 0012 airfoil assumed) 

0.0087 

ii 

7. rotor system Lock number 

2.6 


Table 4.1.1 System Parameters Used in Optimal Control Study 
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speed region is usually described by three points: the high hover point, the low- 
hover point, and the knee or the highest speed point. The low speed region of the 
H-V restriction curve indicates that autorotation from hovering is the most critical 
case. 

In the restricted high speed region there is insufficient clearance between the tail of 
the helicopter and the ground to allow flare to decelerate and control rotor RPM 
prior to ground contact. Since we model the helicopter as a point mass, we cannot 
study the high speed region. The restricted low speed region is of primary interest 
in this work. 

The restricted regions of the H-V diagram are traditionally determined by flight test. 
Using a build-up technique, the pilot performs a series of simulated engine failures 
starting at entry conditions expected to be well outside the restriction regions. 
Subsequent entries are made at more critical conditions until the pilot feels a safe 
landing could not be performed from conditions more critical than the last. These 
tests are performed to establish limits at enough altitude and airspeed combinations 
to allow construction of the H-V diagram. A pilot reaction time of two seconds 
following engine failure is usually employed, along w-ith landing at groundspeeds of 
15 Knots or less. 

Another important factor affecting the H-V restriction is the effect of wind. The 
published restriction areas are shown for zero wind conditions. The presence of a 
headwind when performing an autorotational landing will significantly reduce the 
difficulty of the task. Optimal descents of a helicopter from power loss in hover, 
as well as those with initial forward speed, are given in subsequent sections of this 
chapter for cases with negligible wind effect. 


Another factor that affects the restriction region is the technique used to perform the 







I 
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landings. For the tests reported in reference jl2;. the pilot's primary objective was 
to attain zero rate-of-sink and a safe horizontal speed at touchdown. The technique 
was to level the helicopter 1 or 2 feet off the ground with no rate-of-descent and 
then gradually sink to the ground as RPM continued to decrease. During this phase 
of the maneuver, small amounts of aft cyclic were applied to reduce the horizontal 
velocity. This technique yielded consistent zero rate-of-descent touchdowns while 
the horizontal velocity varied somewhat. 

§4.2 Energy considerations 

Consideration of the helicopter energy state aids in understanding the pilot tech- 
niques required during an autorotational landing, and also the optimal control solu- 
tions given in subsequent sections. After an engine failure, the helicopter has three 
sources of energy: (1) the potential energy of altitude, mgH , (2) the kinetic energy 
of the flight path velocity, jml' 2 , and (3) the rotational energy of the main rotor 
|/Rfl, 2 . Here we have neglected the small contribution of the tail rotor. Therefore 
the total energy (TE) of the helicopter after engine failure is approximately given 
by 


TE = PE + KE + RE , 

= mgH + i mV 2 + i/ R n , 2 . (1) 

Not all the energy given in the above equation is available to the pilot. While 
it is true that we can make “complete” use of the potential and kinetic energy, 
the amount of rotational energy available to the pilot executing an autorotational 
landing is limited. This is due to the fact that the rotor must at all times maintain 
an angular speed above some minimum value Cl mtn . Any drop in rotor speed causes 
an increase in the angle of attack on the rotor’s blades. If not checked, this increase 
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in tli‘* angle of attack will ultimately lead to stall (first on the retreating blades) 
of the rotor and the instability associated wdth it. The minimum rotor speed H min 
defines a bound on the maximum usable rotor energy. The maximum amount of 
energy available to the pilot after engine failure is thus given by 


Maximum Usable Energy = mgH - -mV 2 * -Ir(CI^ - n^, n ) . (2) 

2 2 

The importance of the various energy terms in equation (2) varies with both the 
(H,V) combination at the time of engine failure and the physical configuration of the 
vehicle as represented by the parameters m and Ip. For example, the importance of 
the potential energy term increases with the altitude H at which pow’er failure oc- 
curs. Similarly, the rotational energy term becomes more important either when the 
rotational inertia of the rotor system is increased (such as in the HERS case) or when 
the altitude and airspeed at the time of failure decrease. For a given helicopter model 
(i.e. (m, 1r) = fixed), one can expect the use of different pilot techniques with dif- 
ferent (H,Y) entry conditions. Similarly, under the same entry condition, pilots will 
use somewhat different autorotational landing techniques with different helicopter 
models having different (m,/#) values. 

The time rate of change of the total energy of the helicopter in a typical flight 
situation is given by equation (25) of Section (2.3): 

J t (mgH + ^rnV 2 + ^/flU 2 ) = Ps - (P x + Ppro + Ppara) ■ (3) 


where 


Ps = power supplied to rotor shaft , 
Pi = induced power loss , 

Ppro — profile power loss , 

Ppara = parasite power loss , 


( 4 ) 
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where the various power terms are defined in Section (2.3). 

Equation (3) expresses the energy conservation principle according to which any 
excess power supplied by the engines that is not dissipated by the helicopter is 
stored as internal potential, kinetic or rotational energy. 

In the event of engine failure (i.e. Ps = 0), the total power or energy will decrease 
at a rate which depends on the helicopter configuration as defined by airspeed, main 
rotor thrust, angle of attack and rotor RPM. The pilot's task when entering power- 
off autorotative flight is mainly one of “Energy Management”. He has to control 
the energy losses such that the total amount of energy used from the time of engine 
failure to the final touchdown does not exceed the maximum usable energy (MUE) 
of equation (2). This can be achieved through control of the thrust vector during 
descent in order to land the aircraft at the correct level attitude with small vertical 
and forward speeds. At the same time, during the deceleration phase, the pilot 
must prevent the main rotor from overspeeding which would lead to unacceptable 
blade centrifugal stresses. This task requires considerable pilot skill. 

§4.3 Autorotation Landing Techniques used by Pilots 

For the standard OH-58A helicopter, the autorotation technique depends on the 
flight condition before engine failure. At low altitude and airspeed, below the knee 
of the curve, the pilot technique required for a safe landing is to use increased 
collective to reduce the sink rate as the helicopter approaches 10 to 15 feet above 
the ground. At about 10 feet above the ground, the fuselage is leveled and collective 
is increased as the helicopter settles. At higher airspeeds (still below the airspeed 
at the knee) the collective may be reduced somewhat to regain or maintain rotor 
RPM while the helicopter is decelerated using a cyclic flare. In these conditions, 
the rotor does not enter into a true condition of “autorotation” . 
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For higher altitude entries (above the altitude at the knee of the curve) , the collective 
is reduced and the cyclic moved forward to pitch the nose down in order to gain 
some forward airspeed. The objective is to increase airspeed to that required to 
build rotor speed during a decelerating flare, usually close to the speed defined by 
the knee of the curve. An additional objective of attaining a higher airspeed is 
to establish an operating point at the airspeed corresponding to minimum drag, 
thereby maximizing the amount of time the pilot has to select and maneuver to a 
desired landing area. As the speed increases, altitude rapidly decreases until the 
pilot initiates a flare. The purpose of this flare is to first reduce the rate of sink, 
second, to reduce the forward airspeed before touchdown; and third, to maintain or 
increase rotor RPM. As the helicopter slows down, forward cyclic is applied to level 
the fuselage at the proper landing attitude and, simultaneously, the collective is 
raised to cushion the landing. The techniques described above are taught to Army 
pilots and are described in the OH-58A Flight Manual 59 . 

From an energy point of view, the maneuver of the helicopter, from pilot recognition 
of enginefailure to touchdown, can be divided into three phases; entry, descent, 
and flare. The entry phase consists primarily of arresting the RPM decay of the 
main rotor by lowering the collective pitch angle. This reduces the induced power 
loss P, and established an angle of attack distribution along the blade that results 
in aerodynamic autorotative forces that maintain RPM within the desired range. 
This is especially critical for entry conditions with large induced power loss such 
as when the helicopter is in hover. However, if in hover very near the ground, the 
collective can not be lowered and the rotor RPM will decay immediately. Under 
these conditions, there will not be a steady-state descent phase, and the pilot must 
rely on judicious timing to extract the last available rotor energy to arrest the sink 
rate as ground impact approaches. 
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During the steady-state descent phase, air flows upward through the rotor disk. 
The increase in angle of attack on the rotor blades offsets the drop in the collective 
pitch angle. Total aerodynamic force is increased and inclined forward to establish 
equilibrium. The maneuver in this phase typically consists of a longitudinal cyclic 
“push over'’ to a nose down attitude, followed by acceleration to a desired steady 
descent speed (usually from 50 to 75 knots depending on the helicopter and its 
gross weight). Given enough altitude, steady-state autorotation at this speed can 
be established. Angle of descent is normally 17 to 20 degrees from the horizontal 
(compare to 2 or 3 degrees used in powered descent). The center of attention of the 
pilot in the steady state phase is on maintaining the equilibrium flight through uses 
of the collective and cyclic pitch. In the steady descent phase, some of the potential 
energy of the vehicle has been converted into kinetic energy of the helicopter fuse- 
lage. In effect, the pilot gives up altitude at a controlled rate in return for energy 
to maintain the rotor at a constant RPM. 

During the flare phase, the pilot must reduce airspeed and sink rate just before 
touchdown. Both of these actions can be accomplished by moving the cyclic control 
to the rear. The rearward oriented rotor disk (or TPP) allows a larger volume of 
air to flow through it, resulting in an increase in the total lifting force. The larger, 
rearward pointing thrust will reduce both the airspeed and sink rate. At the same 
time, an in-plane force component is created which causes the rotor to accelerate 
so that some of the lost RPM can be regained. In this phase, kinetic energy of the 
vehicle is converted into rotor energy which, in turn, is used to generate lift. During 
the flare, collective control application may be necessary to prevent rotor overspeed 
or to slow the vertical descent rate even further. Finally, the collective is raised to 
convert the stored rotor energy into lift which further cushions the landing. The 
center of attention of the pilot in the touchdown phase is split between the raising 
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of the collective to cushion the landing and the leveling of the fuselage. The latter 
is performed to avoid striking the tail rotor just prior to touchdown. 

The approximate variations of the various energy terms in the entry, steady descent, 
and deceleration/touchdown phases are illustrated in Figure (4.3.3). 


§4.4 Optimal landing of a Helicopter Initially in Hover 


The optimal descent of a helicopter after power loss in hover was considered by 
Johnson [15]. Instead of having a hard bound on the thrust coefficient Cj- Johnson 
chose to reflect this limitation on the thrust coefficient in another way. The profile 
power loss of the rotor was modified to include a “penalty'” term which increases 
sharply as the loading is raised above the stall limit: 


Vo = " °Cd'l - ( 




) 2 + ( 


kr 

a 

(^F) stall 


) A '*j(l 


4.6/i 2 ) 


(4) 


where ';'(££) , laW is the rotor stall limit and N t is an arbitarily selected large number 
(e.g. JV* =s 20). 


The last term in the square bracket of (4) is the “penalty” term. Its function is to 
increase the profile power loss manyfold w r hen the loading is above the stall limit. 
Its contribution to the profile power loss is insignificant when the rotor operates 
below the stall limit. 


The main advantage of Johnson’s approach is that it avoids the use of path equality 
constraints from the constrained optimization problem of Section (2.9). As a result, 
an ordinary gradient-type, steepest descent algorithm can be used to solve the 
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TE : Total Energy of the Helicopter 

PE : Potential Energy of the Helicopter 

KE : Kinetic Energy of the Helicopter 

RE : Rotational Energy of the Helicopter’* Rotor 

Figure 4.3.3 Variations of Energy Terms in Phases of an 
Autorotational Landing of a Helicopter 


I 



4.4 Optimal landing of a Helicopter Initially in Hover 


125 


resultant TPBVP. The disadvantage of the penalty function approach is that the 
(^F)stali limit is usually exceeded (by as much as 13 percent) over the last 20 percent 
of the trajectory as “ground contact” is approached [15]. 

4.4.1 A Pure Vertical Descent Path 

This section establishes the fact that a pure vertical flight path from power failure in 
hover satisfies the first order necessary conditions. (However, the possible existence 
of a non-vertical flight path with high entry height is suggested by Figure (4.5.27) 
of Section (4.5).) We define a Hamiltonian function H based upon the problem 
formulation of Section (2.0). 

— H = Ai [<?o - motiixl - m 0 fx i Vxj 2 x 2 2 ] 

T f 

+ \ 2 \mQU 2 x\ - mo/x 2 \/zi 2 -t x 2 2 ( 5 ) 

t A3]— tox|(co — Ay^uj 2 — u 2 2 )j — A4 O.IX1! 

Since neither the Hamiltonian function nor the terminal payoff <f> are functions of 
X4 (t?trtical h«ght), we can easily establish the fact that A< is a constant (this is not 
true if ground effect has been included in the formulation; in which case H will be 
a function of height). Knowing this, the Hamiltonian function defined in equation 
(5) can be simplified to 

= Ai|0o - moU!x| - mo/xj V*i 2 + x 2 2 j 

-+ A2[mot/2X3 - mo/X 2 V*l 2 -i- X2 2 ] ( 6 ) 

-f A 3 [— 10 X 3(00 + \/t/i 2 + u 2 2 A)] - 1 - A^O.lxi] • 

ihe optimal controls u\ and 112 are obtained from the optimality conditions of 

h U2 


H Ul — 0 , 


= 0 . 


( 7 ) 
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Using equation (6), these optimality conditions become: 

Aim 0 x 2 - *0^313;-^====== A + A u , \/ui 2 - u 2 2 . = 0 , 

. _ . (si 

A 2 m 0 x 2 -f t 0 A 3 x 2 j — - = — A -r A„ 2 v «i 2 t «2 2 i - 0. 

yu } z - «2 

From the Euler-Lagrange equations, the partial differential equation of A 2 is given 
by 



= -A]m 0 /( - TlJ2 =) - A 2 m 0 / y xj 2 - x 2 2 
V x j 2 - x 2 2 

A 3 i 0 x 2 A I2 \/ui 2 + u 2 2 , . 



( 9 ) 


Using the optimality condition of ( 7 ), and the Euler-Lagrange equations, the optimal 
controls ui and t/ 2 can be solved for simultaneously with the Lagrange multipliers A, 
where i = The analytical solution of the problem cannot be easily obtained. 

However, if we guess that the u 2 component (i.e. the horizontal component) of the 
optimal control U (£) is zero, 

U 2 =0, (10) 

we can proceed to show that the resultant solution of the Lagrange multiplier A 2 
satisfies the optimality condition of (8-2). 

Substitution of (10) into the equation of motion of the normalized forward speed 
x 2 gives 

X2 = r/m 0 /x 2 \/x i 2 + x 2 2 . (11) 

If we integrate this equation forward (in time £ ), using the initial condition of x 2 (0) 
= 0 (i.e. hover condition), the result is 


x 2 (£) = 0, for 0 < £ < 1 . 


(12) 
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Knowing that both uj(£) and X 2 U) are zero, the inflow ratio A defined by equation 
(43) of Section (2.4) can be simplified to 


A — —0.01 kofly‘11] , 

*3 


(13) 


which is neither a function of i <2 nor x 2 . Also, depending on whether it lies within 
the vortex-ring state or the momentum theory state, the induced velocity parameter 
fj is given by equation (14) of Section (2.3). Under assumptions (10) and (12), // 
is now only a function of xj , since ±2 is identically zero: 

. _ Po J i 

' (14) 


*2 


0. 


Using equation (14), we can easily prove that both (//)„ 2 and (//) l2 are zero. 
Therefore, we have from equation (13) that the inflow ratio A is not a function of 
either 1/2 or X 2 : 


K 2 U) = 0, 
W0= 0 . 


(15) 


Making use of the fact that both X 2 U) and A l2 (£) are identically zero in equation 
(9) ghre*.tfie following homogeneous differential equation in A 2 : 

= \T f mofxi('y\ 2 ■ (16) 


Since the terminal payoff <t> is not a function of X 2 (because 12 (f) = 0 for 0 < e< 1 ), 
the terminal value of A 2 is zero. The solution of the homogeneous equation (16) is 
therefore 

A 2 (f) = 0, for 0<f<l. (17) 

Finally, back substitution of equations (10), (15-1), and (17) into (8-2) establish the 
fact that H U2 = 0. Therefore the original assumption of t/ 2 = 0 produces conditions 
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that satisfy the optimality condition of H U2 — 0. Optimal descent of a helicopter 
after power loss in hover as formulated in this section is thus a pure vertical flight 
path. The same conclusion was reached by Johnson in ;15 . 

In practice, a small amount of forward speed is usually required, both to avoid a 
vertical descent into the helicopter's own rotor wake, and also to keep the landing 
point in sight. 

4.4.2 Interpretations of Optimal Control Results 

Both flight tests and analysis indicate that an initial forward velocity greatly im- 
proves the autorotational characteristics of the helicopter. Therefore optimal control 
results are first given for the most critical case of descent from power loss in hover. 
This case was found in Section (4.4.1) to involve only a pure vertical descent. With 
this knowledge, we can eliminate the states 12 (forward speed) and X 5 (horizontal 
distance), and the control U 2 (horizontal component of the thrust coefficient, Cj x ) 
from the problem formulated in Section (2.9). The simplified problem is given in 
Section (3.2) and consists of only three state variables (vertical sink rate, vertical 
height, and angular speed of rotor), and one control variable Cj, (= Ct)- 

The simplified problem was solved here using the Sequential Gradient Restoration 
Algorithm. The Lock number of the OH-58A helicopter considered is 2.61, which 
corresponds to the heaviest blade inertia of 672 slug- ft 2 . The entry heights of the 
helicopter at the time of power loss, given in Figure (4.4.4), vary from 25 to 250 
feet above ground level (AGL). 

Figures (4.4.5) to (4.4.8) present in detail the optimal solution of power-off descent 
from hover from an initial altitude of ho = 100 feet (Case (3) in Figure (4.4.4)). 
The flight time for this case is 4.9 seconds. Figure (4.4.5) gives the collective pitch 
control and the thrust coefficient required as functions of time. Note that the initial 



ENTRY CONDITION 1 2 3 4 5 

ENTRY HEIGHT (ft.) 25 SO 100 200 250 


Figure 4.4.4 Entry Heights of a Helicopter with Power Loss 
in Hover 
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value of the thrust coefficient just before engine failure is given by 


Ct pxRHMt ) 2 ' 

W 

“ 2 ' 

3000 

“ 0.002378 7r (17.65) 2 (37 * 17.65) 2 ' 

= 3.0226 » 10 -3 . (18) 


Therefore the initial values of the inflow ratio A. and the collective pitch angle 
(at 75 percent chord) #75 are given by 


Ct 

— = 0.063 , 
0 

(19) 

[0 1^ 

11 


= 0.03S9. 

(20) 



= 7.12 degrees . 

(21) 


The results given in Figure (4.4.5) agree generally with those from reference '12 . 
The initial drop in the collective is followed by a gradual increase for flare. In Figure 
(4.4.5), the collective flare begins at about one second after engine failure, when the 
helicopter is about 80 feet above the ground. Thereafter, the collective pitch starts 
to increase, and its rate of change reaches a maximum of 6 degrees per second when 
the helicopter is at about the mid-point (in time) of its travel. Because of the hard 
bound on the thrust coefficient at 1 ),<„// = 0.15, the increase in collective pitch 
(or the thrust coefficient) levels off towards the end when touchdown is imminent. 
The thrust coefficient is on the bound during the last one second of travel when 
the helicopter is at about 4 feet above ground level. The results also demonstrate 



Figure 4.4.5 Optimal Time Variations of Thrust Coefficient 
and Collective Pitch Control, ho = 100 feet 
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the effectiveness of the numerical algorithm used in the enforcement of the path 
inequality constraint on the thrust coefficient. 

Figures (4.4.6) and (4.4.7) show the variations of the vertical height, vertical sink 
rate, and the angular speed of the rotor as functions of flight time. 

When the rotor starts to descend, the flows inside and outside the slipstream (cre- 
ated when the helicopter was in the hovering state) are in opposite directions. 
Therefore, from hover to any subsequent windmill brake state, the flow has to first 
go through the vortex-ring state. Figure (4.4.8) shows the variation of the induced 
velocity with the rate of descent in the various states of the rotor in vertical descent. 

The normal hover state and the windmill brake state are well described using mo- 
mentum theory because a definite slipstream exists in these states. By convention 
( see pp. 99 of reference 54 ). the vortex-ring state is defined by P = T(V — u) > 
0 (velocity terms are defined positive in the upward direction), so that the power 
extracted from the airstream is less than the induced power loss. The vortex-ring 
state is also characterized by large recirculation and high turbulence. The region 
with P = I(V -r i/) < Ois called the turbulent wake state. Here, there is a net 
extraction of energy from the airstream. The rotor in this state experiences some 
roughness due to turbulence, but nothing like the high vibration in the vortex-ring 
state. Note that the empirical fairing in Figure (4.4.7) intersects the ideal autoro- 
tation line V ■+■ v = 0 at V/i/^ = 1.6 (see also Figure (3.2.2)). 

In the optimal descent of the helicopter, the rotor is operating in the windmill state 
when the vertical sink rate is increasing (except for a brief period immediately 
after engine failure) . It then operates in the vortex-ring state when the velocity is 
decreasing. Maximum rate of descent is of the order of 2200 fpm and occurs when 
the helicopter is at about 50 feet above ground level. The touchdown sink rate is 


VERT 0E8C«nilL|FT/SEC ) VERTICAL DISTANCE (FT.) 
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Figure 4.4.0 Optimal Time Variations of Vertical Height and 
Vertical Sink Rate [Case (3),h 0 = 100 feet] 
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Figure 4.4.7 Optimal Time Variation of Rotor RPM 
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OF POOR QUALITY 



Figure 4 . 4.8 Variation of Induced Velocity with Rate of Descent [54] 
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of the order of a few feet per second. 

Figure (4.4.8) shows the variation of the rotor RPM with the flight time. The initial 
rate of the RPM decay is about 7 rpm per second. This low rate is due to the drop in 
the collective pitch immediately after power loss. The rate of RPM decay gradually 
increases and reaches a rate of 36 rpm per second at touchdown, when rotational 
energy stored in the rotor is being traded for more lift to cushion the impact. The 
RPM at touchdown is 247, which is about 70 percent of the nominal angular speed. 

4.4.3 Comparison with Flight Data 

To investigate the effect of rotor inertia on autorotational landing performance, a 
flight test program with more than 100 autorotational entries and landings was 
conducted during the HERS flight test programs (cf. (12,57]). Results obtained for 
autorotation landings of HERS from power loss in hover are given in Section (2.4) 
of reference [12:. An optimal flight path was not flown in these tests, of course, 
and in practice, in order to avoid a vertical descent into the rotor’s own wake, some 
forward speed and cyclic flare were involved. Nevertheless, a comparison of flight 
data with results calculated by the optimal programs will shed some light on the 
differences in the piloting techniques involved. 

However, there are difficulties associated with these comparisons. First, it is found 
that the engine torque involved in these flight tests does not decrease to zero im- 
mediately after the throttle is closed to simulate the engine failure. Instead, the 
engine torque pressure decays exponentially with a time constant of approximately 
one second. Second, the grid camera that was used to record the time history of 
the vertical height during these flights covered only up to about 100 feet above the 
ground. The initial portion of the maneuver for a simulated engine failure for a 
hover at 300 feet altitude was therefore out of the range of the camera. We are 
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forced to use results from the only other case, which is a low altitude autorotation 
from an initial height of 50 feet AGL. Third, there is no record of the vertical sink 
rate over time. Finally, for flight tests reported in jl2], the pilot’s primary objective 
was to attain a zero rate of sink at touchdown and to accept any safe horizontal 
velocity. The technique used is to level the aircraft 1 or 2 feet off the ground with 
no rate of descent. The helicopter then gradually sinks to the ground while small 
amounts of aft cyclic are applied to reduce the horizontal velocity. The technique 
consistently yields a zero rate of descent at touchdown but the flight time is usually 
much longer than that calculated from the optimal program. 

Because of the second difficulty mentioned above, we only compare results from a 
simulated engine failure for a hover at 50 feet AGL. The difference in flight times 
obtained from flight data and from calculation make direct comparisons of time 
histories of the rotor speed and the collective pitch control meanningless. The 
problem can be overcome if we instead use height as the independent variable and 
compare plots of rotor speed and collective pitch versus height. Figures (4.4.9) 
and (4.4.10) present the results of these comparisons. Table (4.4.2) summarizes the 
conditions under which the comparisons are made. 

Figure (4 .4 .^compares time histories of the collective pitch calculated by the opti- 
mal program (we shall from here onward call it the computed result) and recorded 
from flight. The comparison is qualitatively good. Two distinct features emerged 
from the comparison. First, relative to the smooth (with an almost uniform increase 
in the collective pitch) computed result, the flight data are wavy. Over two intervals 
in height, first from 5 to 20 and then from 35 to 45 feet, the recorded collective 
pitch remains practically unchanged. The pilot has, perhaps subconsciously divided 
the autorotational landing into phases within each of which his center of attention 
changes. In the interval from 5 to 20 feet, his concentration is on stopping RPM 
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I 



Optimal 

Program 

Flight 
Data [12] 

Cross Weight (lb.) 

3000 

3048 

Rotational Inertia per 
blade (slug- ft 2 ) 

672 

672 

Wind Condition (Knots) 

0 

<3 

Entry Height (ft.) 

50 

50 

Collective Time 
Delay (sec.) 

0 

0 

Flight Time (sec.) 

3 8 

8.1 

Use of Cyclic Pitch’ 

No 

Yes 

Engine Condition 

Engine 

Exponential Decay 

Seizure 

with one second 
time constant 


Table 4.4.2 Conditions Used in Comparison of Calculated 
Results with Flight Data [12] 


























Coiltcttvt Pitch (d*0.) 
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O : Flight Data [12] 


Figure 4.4.0 Comparison of Flight Data with Optimal Results 
: Collective Pitch Control 
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O : Flight Data [12] 


Figure 4.4.10 Comparison of Flight Data with Optimal Results 
: Rotor RPM 
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decay by reducing the collective pitch control. When the helicopter is at about 
10 feet above the ground (interval from 35 to 45 feet), his attention has shifted 
to the use of aft cyclic to reduce both the forward speed and vertical sink rate. 
Preparation is made at this stage for the subsequent raising of collective control to 
cushion the touchdown. It is in such a rapidly changing, often confusing situation 
that an automated control program, which pays attention to “all” aspects of the 
landing maneuver, from the begining till the final touchdown, could be of assistance 
to inexperienced pilots. 

Second, the initial drop in the collective pitch recorded in flight data is only about 
two degrees. This is substantially less than the computed result of five degrees. 
Thereafter the time history of the collective pitch lies consistently below the com- 
puted result. These observations might be explained as follows. In flight tests, the 
engine torque does not drop to zero immediately after engine failure. Instead, it 
decays exponentially with a 1-second time constant to a level which is 10-15 percent 
of the full power output torque. Therefore the residual power available on the rotor 
shaft at the begining of the autorotational descent is significant. This may explain 
why the initial drop in collective pitch is not as much as the computed result. Af- 
ter the initial lowering of the collective pitch, the pilot applies a small amount of 
forward cyclic to pitch the helicopter nose down and attain about 5 knots ground 
speed. This forward motion reduces the induced power consumption of the heli- 
copter which explains why recorded data lie consistently below computed results. 
The gain obtainable with the use of forward cyclic is of course not without a price. 
Aft cyclic must be applied in the landing flare phase of the descent to reduce the 
forward speed of the helicopter to zero ground speed at touchdown. In addition 
to the increased induced power loss associated with the use of aft cyclic, the total 
flight time of the autorotation descent is lengthened. 
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Figure (4.4.10) compares the recorded time history of rotor speed with that calcu- 
lated by the optimal program. While one should not look for too much correlation 
here, the comparison is surprisingly good. The rotor speed at touchdown is 220 
rpm. This value is somewhat lower than the computed result of 268 rpm (see also 
Table (4.4.3)). The total amount of energy spent in the optimal program is 

Energy Used (Optimal Program) = mgH -t- -//?(fl? - fly) 

= 3000 x 50 - -(672 x 2)(37 2 - 28 2 ) ( 22 ) 

2 

= 54.3 x 10 4 ft. - lb.wt. 

This amount is about 24 percent less than that used by the pilot in achieving the 
same zero rate of sink at touchdown from a 50 feet hovering throttle chop: 

Energy Used (Pilot) = 3000 x 50 + ^(672 x 2)(37 2 - 23 2 ) , 

2 (23) 

= 71.4 x 10 4 ft. — lb.wt. 


4.4.4 Effects of Entry Height 

Figure (4.4.11) compares the time history of the collective pitch control calculated 
by the optimal program at three different entry heights of 50, 100, and 200 feet. 
Since the computed optimal flight time varies with the entry height, we change the 
independent variable from time to the normalized time (£), where £ is defined as 

e = p, (24) 

*/ 

where tf is the flight time. Computed time histories of the thrust coefficient are 
compared in Figure (4.4.12) 

Figure (4.4.11) shows that the optimal collective pitch control program changes with 
the entry height. From the 7.2-degree collective pitch used before engine failure, 
the initial collective pitch reduction is in the range of 5-6 degrees for all the entry 


Flight 
Data [12] 


Optimal 

Program 


Time of Touchdown (sec.) 


Max. Rate of Descent (fps. 


Touchdown vertical speed 
(fps.) 


Rotor Speed at Touchdown 
(rpm) 


Minimum Collective Pitch 
(degrees) . . 


Col lective Pitch at 
Touchdown (degrees) 


8.1 


3.8 


IS 


20 


221 


268 


2.0 


5.2 


14.2 


14.8 


Table 4.4.3 Comparison of Optimal Results with 
Flight Data, ho = 50 feet [12] 
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0.25 


0.50 


0.75 


100 


Normalized time 


ENTRY CONDITION 

2 3 4 

ENTSY HEIGHT (ft . ) 

50 100 200 

FLIGHT TIME (sec.) 

3.8 4.9 6.1 


Figure 4.4.11 Time Variations of Collective Pitch Control 
from Different Entry Heights 
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ENTRY CONDITION 

12 3 

ENTRY HEIGHT (ft.) 

25 50 100 

FLIGHT TIME (sec.) 

ma 


Figure 4.4.12 Time Variation of Thrust Coefficient 
from Different Entry Heights 
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heights considered. The collective is lowered to minimize RPM decay. At the lowest 
entry height of 50 feet, the collective is immediately raised again after the initial 
drop in collective. At the higher entry heights of 100 and 200 feet, the lowered 
collective pitch is maintained for some time before it is raised to reduce the rate 
of descent. This delay allows a build-up in the rate of descent of the helicopter, 
but with the high rotor inertia, subsequent collective increase provides a sufficient 
deceleration to allow a safe landing. 

4.4.5 Most Critical Entry Height 

Figure (4.4.12) shows the time history of the thrust coefficient calculated by the 
optimal program at three different entry heights of 25, 50, and 100 feet. The curve 
for case (4) (w’ith entry height of 200 feet) lies slightly to the right of that for case 
(3) and has been omitted in the interest of clarity. 

As the entry height is increased, Figure (4.4.12) shows a corresponding increase 
in the amount of (normalized) time that the thrust coefficient spends on its stall 
limit. This trend continues till the entry height of 100 feet is reached. Beyond the 
“critical” height of 100 feet, any further increase in entry height reduces the time 
that the thrust coefficient is on its stall limit. Since prolonged usage of the maximum 
available thrust is an indication of an irrevocable extraction of the rotational energy 
from the rotor system during the autorotational landing maneuver, the observed 
results suggest the existence of a “Most Critical Entry Height” (MCEH). Above 
or below- this critical height, which in the present case is around 100 feet, the 
autorotational landing procedure becomes progressively easier to execute. 

Further confirmation of this hypothesis can be found in the time history of the 
rotor angular speed in Figure (4.4.13). Since the rates of descent at touchdown for 
all entry heights considered are practically zero, we cannot use them to compare 
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the criticality of the landing". However, we have learned from Section (4.4.3) that 
the lower the rotor RPM at touchdown, the larger is the amount of energy (of the 
helicopter) used in the landing maneuver. One easy way to show the existence of 
the MCEH is to plot the rotor speed at touchdown against the corresponding entry 
height of the autorotational descent. Such a plot is given in Figure (4.4.14). 

Figure (4.4.14) shows that the MCEH (denoted by h cr ) is at around 100 feet where 
the rotor speed at touchdown is a minimum. Therefore the “Most Critical Entry 
Condition" is one with engine failure in hover at h CT . This conclusion agrees with 
flight test results recorded in reference 'll . 

§4.5 Optima) Landing of a Helicopter in Forward Flight 

The autorotational landing of a helicopter in forward flight has been formulated as 
an optimal control problem in Section (3.3). Unlike the case when power loss occurs 
in hover, the present case involves more than a pure vertical descent. As such, we 
have to consider the “full state” problem w'ith five states 

Xi = normalized vertical velocity , - 

xj = normalized horizontal speed. 

X3 = normalized rotor speed , (24a) 

X 4 = normalized vertical height , 
x$ = normalized horizontal distance . 


and two controls 

u i = vertical component of thrust , 

«2 = horizontal component of thrust . 


(25) 


If the helicopter is not constrained to land at a particular spot on the ground, the 
6tate x& may be removed from the list in (24). Optimization is then performed on a 




Rotor spe«d (rpa) 
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Normalized time 




Time Variations of Rotor Speed at Different 
Entry Heights 



impact Velocity ppsj 
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problem with only four state? and two controls. After the optimal solution has been 
computed, information on the horizontal distance (x 5 ) may be recovered through 
the forward integration of the kinematical relation: 

X 5 (C) = f l \0.lT f x 2 U) d'. (26) 

Jo 

where £ is the normalized flight time. 

4.5.1 Flight Program 

Once again, the High Energy Rotor System (HERS) described in Section (4.1) is 
used as the basis for the analytical model. However, we note in reference [12] that 
the highest inertia rotor (with -y = 2.61) allows complete elimination of the height- 
velocity restriction curves. For the mid- inertia rotor (with -y — 3.19), only a small 
restricted region remains between 75 to 125 feet altitude with airspeed less than 5 
knots. Therefore it is more interesting to investigate optimal landings of the HERS 
with its lowest inertia rotor. We have chosen, rather arbitarily, a rotor with Lock 
Number *) of 4.38 which corresponds to a rotor inertia (per blade) of 400 slug-ft 2 . 
The low-speed height-velocity restrictions for the HERS with rotor inertia at three 
different levels are given in Figure (4.5.15) along with the diagram for the standard 
OH-58A helicopter. 

Entry conditions studied are summarized in Figure (4.5.16). We have chosen to 
analyse entry conditions that are within or close to the height-velocity restriction 
curve of the HERS with its low-inertia rotor. Note also that entry conditions 1, 2 and 
6 are at the same altitude of 100 feet but with different forward speeds. Similarly, 
entry conditions 1, 5 and 4 have (approximately) the same forward speed of 8 knots 
but are at different entry heights. In this way, effects of both the forward speed and 
altitude at entry can be studied. Optimal landings of the HERS studied at these 



Figure 4.5. IS Comparison of HERS Test Results with Standard OH-58A [12] 
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Figure 4.5.16 (H-V) Entry Conditions 
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entry conditions are without any terminal constraint on the horizontal distance. 
The only exception is entry condition 5, where the helicopter is constrained to land 
at various horizontal distances from the point of engine failure. 

4.5.2 Interpretation of Results 

Figures (4.5.17) to (4.5.19) present in detail the optimal solution of power-off descent 
from an altitude of 100 feet and at 12-knots forward speed (entry condition 1 in 
Figure (4.5.16)). Figure (4.5.17) gives the collective pitch control and the thrust 
coefficient as functions of time. Since the parasite drag on the fuselage of the 
helicopter in level flight is small compared with the weight of the vehicle, equation 
(18) is still approximately true and the value of the thrust coefficient divided by 
the rotor solidity before engine failure is given (approximately) by 0.063. 

The induced velocity in level flight can be found from the momentum quart ic: 

v 4 - 21’ £> 3 sin q -h V 2 P 2 - (^) 2 = 0. (27) 


where the normalized induced velocity (i>) and the normalized flight path velocity 

(V) ift^ivenby 



(28) 


In level flight, the angle which the flight path velocity V makes with the TPP, a, 
is approximately zero. Equation (27) can thus be simplified and the value of £> is 
given by the following equation 
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The value of the induced velocity at hover. Uh- is given by 


- — 

Uh \ 2pA ' 

= 25.38 fps. 


(30) 


Values of I ' and v can then be computed from equations (28-2) and (29) respectively. 
Their values are 

V = 0.798. 


(31) 


v = 0.855. 

Therefore the inflow and advance ratios are given by 

p = 0.031 , 

(32) 

A = 0.0332 . 


The collective pitch used before engine failure can then be computed from the 
following relation 

_ (i + f^ 2 )( 6 -£) + |Mi-^ 2 ) 

73 (l_ M 2 + ^ 4) 

. 6Cj- 3 (33) 

= -T- — A . 

aa 2 

= 6.63 degrees. 

As expected, this value is smaller than the value used in the hover case (7.72 de- 
grees). 

The optima] solution is similar to that found in the hover case. Initially the collec- 
tive is lowered to maintain rotor speed, and nose-down cyclic is used to maintain 
airspeed. The results of these control actions are an initial build-up of the verti- 
cal sink-rate as well as an increase in the forward airspeed. However, the rate of 
RPM decay is being curtailed. In fact, the angular speed of the rotor stabilizes at 
around 332 rpm until a rearward, nose-up cyclic flare is used to slow both the rate 
of descent and forward airspeed. 


COLLECTIVE WCH (D€G| 
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Figure 4.5.17 Optimal Time Variations of Thrust Coefficient 
and Collective Pitch Control 
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Figure 4.5.18 Optimal Time Variations of Horizontal and 
Vertical Velocity 
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Figure 4.5.10 Optimal Time Variations of Rotor RPM and 
Helicopter's Altitude 
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Because of the bound on thrust coefficient at (^-) — 0.15, the steady increase in 
the collective pitch levels off towards the end when touchdown is imminent. The 
thrust coefficient stays on its bound during the last second of the travel when the 
helicopter is very near to the ground. Once again, the path inequality constraint on 
the thrust coefficient has been effectively enforced by the optimization algorithm. 

Nose-up cyclic is initiated 2 seconds after the engine failure when the helicopter is 
at about 50 feet above the ground. This occurs later than was the case for hover, 
yet there will be sufficient kinetic energy available for the landing flare on account 
of the initial forward speed. With a rearward tilt of the TPP, the projected area 
of the TPP in the direction perpendicular to the air flow has increased. More air 
now flows through the rotor disk. The resultant increase in the angle of attack of 
the rotor blades increases the thrust. This effect, together with the steady increase 
in the collective pitch, bring about the desired reductions in the forward speed and 
rate of descent approaching touchdown. 

The time variation of rotor speed is given in Figure (4.5.19). The figure clearly shows 
that some of the lost RPM can be regained if a nose-up cyclic flare is initiated with 
sufficiently high helicopter forward speed. Effects produced by a rearward (nose- 
up)cvclic flare are illustrated in Figure (4.5.20). Because of the increased lift force 
(i.e. L >> D in Figure (4.5.20)), the resultant thrust vector R is tilted forward 
and produces a net force component in the plane of the rotor’s blades. This force 
component generates a torque that accelerates the rotor. During the touchdown 
phase, the collective is raised to cushion the impact. The rotor speed is reduced to 
225 rpm at touchdown, indicative of extracting most of the available rotor energy. 
The touchdown is made at a near zero horizontal speed. 


4.5.3 Comparison with Flight Data 
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In this sub-section, techniques used by pilots are compared with those computed 
from the optimal program. We use flight data obtained from a simulated engine 
failure of the low-inertia version of the HERS at the “knee’* of the height velocity 
restriction curve ( 12 . Flight data obtained from this entry condition (with 115 feet 
altitude and 45 knots forward speed) are compared with those computed by the 
optimal program for the entry condition at 100 feet altitude and 38 knots forward 
speed. Therefore, one must note the small differences in both the entry height and 
speed when making comparisons. Additionally, because of the strong emphasis on a 
zero vertical sink-rate at touchdown, pilots tend to neglect the horizontal component 
of the velocity in their landing maneuver (both military and civil criteria permit 
this). In this particular case, the touchdown horizontal speed of the helicopter is 
20 fps. This pilot priority is quite different from the cost function used in the 
formulation of the optimal control problem. Other differences in flight conditions 
are summarized in Table (4.5.4). 

Because of the difference in flight times obtained from flight data and that from 
calculation, direct comparisons of the time histories of the collective pitch and the 
rotor angular speed are not meaningful. The previously used technique of converting 
the independent variable from time to height is not applicable here because of the 
difference in entry heights. Therefore we choose to use the normalized time (£, see 
equation (24)) as the independent variable in making comparisons. Figures (4.5.21), 
(4.5.22) and (4.5.23) present results of these comparisons. 

Figure (4.5.21) compares time histories of the collective pitch calculated by the 
optimal program and from flight data. The comparison is qualitatively good. The 
most obvious differences in these time histories are the undulating nature of the 
flight data, as well as the fact that the collective pitch control used in flight tests 
is lower than that calculated. One of the reasons for the observed difference in 
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Optimal 

Program 

Flight 
Data [12] 

Cross Weight (lb.) 

3000 

3040 

Rotational Inertia per 
Blade (slug- ft 2 ) 

400 

323 

Wind Condition (Knots) 

0 

<3 

Entry Height (ft.) 

100 

115 

Entry Speed (Knots) 

38 

45 

Collective Time Delay 
(sec.) 

0 

0 

Touchdown Vertical 
Speed ( fps . ) 

0 

0 

Touchdown Horizontal 
Speed ( fps . ) 

0 

20 

Rotor Speed at 

Touchdown (rpm) 

224 

280 

Flight Time (sec.) 

8.7 

10.8 

Engine Condition 

Engine 

Seizure 

Exponential 
Decay with 
1- second time 
constant 


Table 4.5.4 Flight Test Conditions 


























Collective Pitch (degrees) 
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O Flight Data [12] 


Figure 4.5.21 Comparison of Flight Data with Optimal Results 
: Collective Pitch Control 
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collective pitch controls is the availability of the residual power to the helicopter 
in flight tests. The difference is also partly due to the fact that the touchdown 
horizontal speed of the helicopter in the flight test is 20 fps (as compared with the 
near zero value in the optimal result). With reference to Figure (4.5.18), we can 
see that the calculated vertical sink-rate of the helicopter over the last 1/3 of the 
flight time is practically zero. Over the same period of time, the forward speed 
of the helicopter is being steadily reduced to its near zero value at touchdown. 
Therefore the touchdown condition of the flight test (zero rate of descent and 20 
fps forward speed) has been achieved in the optimal program at around 2/3 of the 
flight time. If we now compare the value of the computed collective pitch at around 
2/3 of the flight time with that from flight data at touchdown, their values become 
comparable. Unlike the flight test, the optimal program will not land the helicopter 
with high forward speed. Instead it will try to reduce the forward speed to as near 
to zero as possible. This is done over the last 1/3 of the flight time where both the 
collective pitch and rear cyclic are used to further decelerate the helicopter. The 
above reasoning helps to explain the difference in collective pitch controls shown in 
Figure (4.5.18). 


Figure (4.5.22) compares the recorded time history of the rotor speed with that 
calculated by the optimal program. The comparison is again qualitatively good. 
Notice that the pilot had placed much greater emphasis on maintaining a constant 
rotor RPM than that reflected in the optimal results. Consequently, the collective 
pitch is maintained at a minimum for as long as possible. This helps to explain the 
discrepancy found in Figure (4.5.22). Since the optimal program decelerates the 
forward speed of the vehicle to near zero while flight test to 20 fps, the computed 
rotor RPM at touchdown (about 230 t pm ) is lower than that found m the actual 
flight test (280 rpm). 



O 0.1 0.2 0.3 0.4 0.0 O a O.T O O 0.0 1 

Normalized Tima 


O Flight Data [12] 


Figure 4.5.22 Comparison of Flight Data with Optimal Results 
: Rotor RPM 
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A comparison of the optimal trajectory with that found in flight tests is shown in 
Figure (4.5.23). Flight data show a 535 ft horizontal distance at touchdow’n, which 
is larger than that from the optimal program (400 ft). The extra distance was used 
by the pilot in the reduction of forward speed through the use of cyclic control. 
In order to maintain the helicopter at its level attitude for the final touchdown, 
only small amount of rearward cyclic was used in the deceleration of the vehicle. 
This pilot actions results in a prolonged touchdown and a longer horizontal distance 
travelled. 


Landing Techniques 


The most critical phase of the autorotation descent of a helicopter is the final 
approach to a landing flare. Figures (4.5.24) and (4.5.25) show the pilot control 
motions and flight profiles from the flight data, and those obtained by the optimal 
program, respectively. Both records are for an autorotation entry condition of about 
40 to 45 knots airspeed, and level flight at around 100 feet altitude. Superimposed 
on the flight profile are reference lines indicating the orientation of the helicopter's 
pitch attitude at one-second intervals. Close spacing of the reference lines indicates 
a slow maneuver while larger spacing indicates a faster maneuver. As the helicopter 
is modelled as a point mass, pitch attitude of the vehicle is not available from the 
optimal program. In Figure (4.5.25), we have substituted for it the orientation of 
the TPP in space. Also, since the longitudinal cyclic cannot be determined from 
the optimal program, it is not shown in Figure (4.5.25). 


As can be seen from the flight data, a cyclic flare is initiated just below 100 feet 
and held to about 20 feet above the ground. This results in a nose-up attitude 


of IS degrees. The collective ts increased slightly during the flare to prevent rotor 
overspeed. Only after reaching about 15 feet does the final collective pull begin. In 
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O Flight Data [12] 


Figure 4.5.23 Comparison of Flight Data with Optimal Results 
Flight Trajectory 
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contrast, the computed results given in Figure (4.5.25) show that the collective is 
pulled earlier at an altitude of about 30 feet and is raised almost continuously until 
touchdown. Rear cyclic is applied at about 70 feet altitude to reduce the horizontal 
velocity. The maximum deviation of the TPP from horizontal during the descent 
is 24 degrees (in the nose-up direction). This angle is reduced progressively as the 
end-time is approached. The return of the thrust vector to the vertical is done 
because most of the forward speed of the vehicle has already been reduced as the 
end-time is approached. It also provides additional lift to cushion the touchdown. 
In an actual landing, the vehicle must be rotated to its landing attitude so as to 
provide sufficient clearance between the tail-rotor and the ground. This particular 
constraint is not included here because of our use of a point mass model. 

4.5.4 Effects of Entry Height 

Time histories of optimal thrust coefficient for entry conditions with approximately 
the same forward speed (8-12 knots), but at three different entry heights of 100, 
230. and 420 feet, are compared in Figure (4.5.26). Since the computed flight time 
varies with the entry height, we change the independent variable from time t to the 
normalized time £. Figure (4.5.26) shows that the optimal thrust program changes 
with the entry height. In all three cases considered, drops from its initial value 
of 0.063 (used before engine failure) to a lower value in order to reduce the initial 
RPM decay. This initial drop in ^ is most severe at the lower entry height of 100 
feet. The reduction is relatively mild at higher entry heights of 230 and 420 feet. 

After the initial reduction, the thrust coefficient then increases monotonically. The 
rate at which is increased is again a function of the entry height. At the lower 
height of 100 feet, there is a relatively rapid increase in the thrust coefficient. Be- 
cause of the upper bound on the rate of increase levels off when the end-time 
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Figure 4.5.24 Flight Data for Autorotation Landing [12] 
[ Entry Condition : 125 feet, 45 knots ] 
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Figure 4.6.25 Computed Optimal Program for Autorotation Landing 
( Entry Condition : 100 feet, 38 knots ] 
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ENTRY CONDITION 

15 4 

ENTRY HEIGHT (ft.) 

100 230 420 

FLIGHT TIME (ssc.) 

5.5 9.3 11.3 


Figure 4.5.26 Variations of Thrust Coefficients with Entry Height 
[ Entry Speed : 8-12 knots ] 
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is approached. The thrust coefficient stays on its stall limit = 0.15) over the 
last second of its travel in the optimal control program. 

The same trend is followed in cases when engine failure occurs at the higher altitudes 
of 230 and 420 feet. However, the larger potential energy which is now available to 
the helicopter makes it unnecessary to have a large increase in thrust coefficient. 
This is because some of the potential energy of the helicopter can be traded for 
kinetic energy, and that the same increase in collective pitch produces a larger 
thrust at a higher forward speed than at a lower speed. Also, throughout the 
complete descent process, the stall limit is never reached. These observations agree 
with those found in the hover cases. In these cases, a “most critical entry height’" 
(MCEH) exists from which the autoroational descent is more difficult to execute 
than from any other entry heights. This concept of MCEH can be extended to entry 
conditions with the same entry speed but different entry heights. As is indicated 
in Figure (4.5.26), with approximately the same entry speed of 8-12 knots, the 
autorot at ion al landing from the entry height of 100 feet is more difficult to execute 
than from any other entry height. 

Further insights from the optimal programs can be gained if we resolve the thrust 
coefficient into its horizontal and vertical components and plot these thrust compo- 
nents against the normalized time. Such a plot is given in Figure (4.5.27). Variation 
of Cr, with the normalized time is qualitatively similar to that of the thrust co- 
efficient shown in Figure (4.5.26). In all three cases considered, the horizontal 
component of thrust coefficient is positive for the first half (in time) of the de- 
scent, indicative of the use of forward cyclic to accelerate the helicopter to a higher 
speed. This forward acceleration, of course, cannot be continued for too long, and 
at around mid-time, rear cyclic is applied to decelerate the helicopter and prepare 
it for the final touchdown. Time histories of the forward speeds obtained with such 
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an acceleration /deceleration optimal control program are shown in Figure (4.5.28). 

In Figure (4.5.28), the forward speed first increases, reaches a peak at around mid- 
time and then decreases to zero at touchdown: the higher the entry height, the 
larger the peak speed. The higher forward speed in turn increases the effectiveness 
of the rearward (nose-up) cyclic flare (i.e. with negative Cj- X ). This explains "hj 
touchdown is made with progressively lower ^ value in Figure (4.5.26) as the entry 
height is increased. 

4.5.5 Effects of Entry Speed 

Time histories of optimal thrust coefficient with entry speeds of 12, 38, and 57 
knots, all at an entry height of 100 feet, are compared in Figure (4.5.29). Entry 
height of 100 feet is selected because of the “criticality" of this height discussed 
previously. Optimal time variations of the horizontal and vertical components of 
the thrust vector are given in Figure (4.5.30). Time variations of forward speed for 
these entry conditions are given in Figure (4.5.31). 

Figure (4.5.29) shows the effect that entry speed has on the thrust program. The 
time variation of in entry condition 1, which is the case considered in Section 
(4.5.4), consists of an initial drop in followed immediately by a steady increase 
that levels off as the end-time is approached. Variations in cases with higher entry 
speeds are qualitatively similar to that of the first case. The higher the entry speed, 
the larger is the kinetic energy available to the helicopter in its execution of the 
autorotation landing. Since the effectiveness of the rearward cyclic flare increases 
with flight path velocity, there is no need for a large increase in collective at the 
end of the maneuver for entries with high forward speed. Hover, a condition with 
zero forward speed, is therefore the most critical case to execute an autorotation 
landing for any given entry height. 


Thrust coefficient 



Normalized time 



Figure 4.S.27 Variations of Horizontal and Vertical Thrust Components 
with Entry Height [ Entry Speed: 8-12 knots ] 
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Figure 4.5.28 Optimal Time Histories of Forward Speed at different 
Entry Heights [ Entry Speed : 8*12 knots ] 
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Normalized time 


ENTRY CONDITION 


ENTRY SPEED (Knots) 

12 38 57 

FLIGHT TIME (sec.) 

5.5 8.7 10.4 


Figure 4.5.20 Variations of Thrust Coefficients with Entry Speed 
[ at Entry Height of 100 feet ] 
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Figure 4.5.30 Variations of Horizontal and Vertical Thrust Components 
with Entry Speed [ at Entry Height of 100 feet ] 
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An interesting feature of the optimal program is observed in Figure (4.5.30). With 
the same entry height of 100 feet, we note the change in the Cj x -program with 
entry speed. If the forward speed at the time of power loss is relatively low (e.g. 12 
knots), forward cyclic must be used to accelerate the helicopter to a higher forward 
speed before the aft cyclic is initiated to slow it down for the touchdown. At an 
intermediate speed of 38 knots, there is but a brief period of forward acceleration. 
At the highest forward speed considered (57 knots), aft cyclic must be applied 
immediately after engine failure in order to decelerate the helicopter in time for the 
final touchdown. Time variations of forward speeds in these three cases considered 
are given in Figure (4.5.31). 

4.5.6 Best Endurance Speed 

The need to accelerate the helicopter to a higher forward speed when engine failure 
occurs at low entry speeds and to decelerate it when power is lost at high speeds 
indicates the existence of some intermediate speed at which the helicopter is neither 
accelerated nor decelerated. To reduce workload, a pilot might prefer to initiate an 
autorotation landing from this speed. 


This “preferred” speed can be better understood if we study the effect of forward 
speed on the energy consumption of a helicopter in level, powered flight. Total 
power consumption is the sum of induced power loss of the rotor, blade profile 
power loss, fuselage parasite power loss, and tail rotor power. Tail rotor power and 
blade profile power are not strong functions of the forward velocity. In contrast, 
induced power loss decreases drastically with increased forward velocity while the 
parasite power loss increases with the third power of the forward velocity. Figure 
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the characteristic “bucket” shape. The speed for best endurance is the speed at 
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FLIGHT TIME (soe.) 

5.5 8.7 10.4 


Figure 4.5.31 Variation of Forward Speed with Entry Speed 
at Entry Height of 100 feet 
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which the power-required curve is at its minimum. Tlu speed for best range, which 
is usually slightly higher than the speed for best endurance, is the speed at which 
the power-required curve is tangent to a line drawn through the origin. This is also 
the point at which the ratio of speed to power (and therefore, of distance to fuel) 
is greatest. 

In autorotation, steady flight can be achieved only by descending. The rate of 
steady autorotative descent is once again a strong function of the forward speed 
of the helicopter. Appendix (A) presents a way to compute the rate of steady 
descent at different forward speeds. The rate of steady descent in autorotation of 
the OH-58A helicopter considered here is shown in Figure (4.5.33). Again, we see 
the characteristic “bucket" shape, indicative of some intermediate speed at which 
the rate of descent is a minimum. The forward speed w r ith the minimum rate of 
descent is the same as the speed for best endurance shown in Figure (4.5.32). This 
speed for the OH-58A helicopter is at around 45-50 knots. 

From the point of view of energy mangement, the optimal control program tries to 
fly the helicopter at a speed which minimizes the power loss. The optimal program 
achieves this by accelerating or decelerating the helicopter to a higher or lower 
forward speed as the case may be. However, we note that results shown in both 
Figures (4.5.32) and (4.5.33) are obtained with the helicopter in a steady state 
condition. Since the optimal results are in the transient state, we cannot compare 
these optimal results with those shown in the figures directly. Nevertheless, it is 
still not too surprising to see that the optimal program has avoided both high and 
low forward speed with very high power consumption. 

§4.6 Optimal Descent of a Helicopter with a bound on Vertical Sinkrate 
In the last section, we presented and interpreted some of the results obtained from 
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Figure 4.5.32 Variation of Power with Forward Speed 
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Figure 4.5.33 Variation of Steady Rate of Descent with Forward Speed 
of an OH-58A Helicopter 
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Figure 4.5.34 Time Variations of Rotor Speed 
[ Entry Speed : 8*12 knots ] 
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Figure 4.5.SS Time Variations of Rotor Speed 
[ Entry Height : 100 feet ] 
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computations of optimal descents of a helicopter in forward flight. The optimal 
descents of the HERS from entry conditions either with a high forward speed or at 
a high altitude result in peak rates of descent that are unacceptably high to pilot, 
and to overspeeding of the rotor during the rearward cyclic flare. 

To remedy this problem, we add an upper bound on the vertical sink-rate of the 
helicopter to the optimal landing problem. The new problem, with an additional 
state variable inequality constraint is given in Section (3.3). An entry condition 
with a 7.7-knot initial speed at an initial altitude of 423 feet was selected for the 
present study. This entry condition was used because the optimal descent of the 
HERS with this particular entry condition results in unacceptably high peak descent 
velocity. Results obtained with the optimal descent of the HERS from this entry 
condition and subjected to a descent velocity bound are presented in this section. 

4.6.1 Optimal Descent of a Helicopter with a Descent Velocity Bound 

If the optimal rate of descent of the helicopter in the unbounded solution exceeds 
the upper bound, it is reasonable to assume that the optimal time variation of the 
descent velocity in the solution to the new problem will resemble that shown in 
Figure (4.6.36). The helicopter starts initially with a zero rate of descent. This rate 
then increases, touches and stays on its upper bound for a period of time. Finally, 
the rate of descent is reduced so as to allow for a soft touchdown. Typically, this 
time variation of the descent velocity consists of three phases. These phases of 
entry, steady descent, and landing flare are shown in Figure (4.6.36). 

The amount of time that the helicopter spends in each of these phases depends 
on the magnitude of the descent velocity bound, [Vp] m0I . Generally speaking, the 
higher the descent velocity bound, the shorter the time spent in the steady descent 
phase. As illustrated in Figure (4.6.37), the steady descent phase ceases to exist 
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Time (seconds) 


(1) Entry Phase : —2= positive, 

[ phase ends with V D = (V D ) max ] 

d Vn 

(2) Steady Descent Phase : "dt = 2ero ' 

c V D = (VdW ] 

dVn 

(3) Landing Flare Phase ; "dt * ne 9 a t:ive, 

[ phase ends with V D s O ] 

Figure 4.6.36 Time Rate of Change of Descent Velocity 
with a [V D ] mat Bound 
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(1) High : The variation is essentially that 

without a descent velocity bound. 

(2) Medium (Vp] max ‘The variation consists of the entry, 

steady descent and landing flare phases 

(3) Low [Vp) m „ : No feasible solution exists in this 

case. 


Figure 4.6.37 Time Variation of Descent Velocity with 
Different Values of [Vp] ma , 




Rat* of Descent 


4.6 Optimal Descent of a Helicopter with a bound on Vertical Sinkrate 


187 


ORIGINAL PAGE IS 

OF POOR QUALfTY 


2400 


20001 


1400 


1200 


800 



“1 

:ons' 

IMSt 

Rotor K 

— 

m • 

^354 

r~ 





• 











ft 



A 

T 



. J. 

^ ... . 

0 



J 

V 





r ™' l 

&©< 

Mm 


l ^ 

V 

0 







oo 







i 

0 

Di 

it* < 

3f R 4 

• fen 

inca 

5 




■ o 

fans 

— Fairing of Reference 

5 

' 



20 40 60 

Airspeed , kn 


80 


100 


Figure 4 .0.38 Variation of Steady Autorotative Sink-rate of the HERS 
with Forward Airspeed 
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when the 'Vy> mQI -bound is chosen too high. In this c<:»e, the variation of the 
descent velocity becomes essentially that without the descent velocity bound. The 
other extreme is when the !tz?j mQ2 'b° un d selected is too low, e.g., at around 1200 
fpm. As can be seen from Figure (4.6.38), the horizontal line .} D] maz = 1200 fpm 
does not intersect the curve which shows the variation of the steady autorotative 
sink-rate of the HERS with its forward airspeed. Therefore, no feasible control 
exists that will allow the helicopter to autorotate at a steady sink-rate of 1200 
fpm. 

An intermediate value of 1800 fpm was selected as an upper bound on the vertical 
sink-rate of the helicopter in autorotative descent. From here on we shall use the 
term “nominal case* to represent the autorotative descent of the HERS from an 
entry condition with a 7.7-knots forward airspeed, at an initial altitude of 423 feet, 
and having a 1800 fpm bound on its vertical sink-rate. 

The time variations of the Cy,, Cy a , and the collective pitch for the nominal 
case are given in Figure (4.6.39)-(4.6.40). With reference to the Cj t - curve, one can 
see a natural division of the autorotative landing maneuver into the three phases 
mentioned earlier. Similar divisions can also be observed in both the time histories 
of the collective pitch and the thrust coefficient. 

An approximate Solution of the Cy, 

During the steady descent phase, when the sink-rate of the helicopter stays close 
to the value [Vp] mor , the vertical component of the thrust coefficient Cy,, can be 
determined using the following approximations. During the steady descent phase, 
the time rate of change of the helicopter’s sink-rate is near zero. By equating the 
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Figure 4.6.30 Time Variations of the Thrust Coefficient and its 
Horizontal and Vertical Components 
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Figure 4 . 6.40 Time Variation of the Collective Pitch 
with a Descent Velocity Bound 
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LHS of the first equation of Appendix (A) to zero, we obtain 


0 = g 0 - m 0 (i/]X§ + fx\ \J x] t- x ]) , 
£ = u,x^/x,\/xf-x2. 


(34) 


Since ^ = 3.02, / = 0.0306, and magnitudes of both xj, X2 and X3 are of the order 
of 1. the second term in equation (34) (i.e., the parasite-drag term, fx\yj~x\ - x]) 
may be neglected. Therefore, we obtain the approximate relation: 

— = 3.02. 

m o (35) 

= UjX^. 


which can be used to estimate the value of the normalized vertical thrust coefficient 
«i- 


If the angular speed of the main rotor does not vary too much during the steady 
descent phase, the normalized rotor speed X3 can be approximated by its mean 
value |x3] mean . Therefore the approximate value of U] in the steady descent phase 
is given by: 


_ go 

^0 (-^3) m eon 
3.02 

“ (0.95) 2 ’ 

= 3.33 . (36) 


where we have used a value of 0.95 for the average normalized rotor speed ((x3) mea n) 
in the steady descent phase (see also Figure (4.6.43)). The computed value of 3.33 
for the vertical component of the thrust coefficient in the steady descent phase 
agrees well with the value 3.7 found in the optimal solution (see Figure (4.6.41)). 

4.6.2 Comparisons between Results Obtained With and Without 
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Figure 4.6.41 Time Variations of the Thrust Coefficients With and 
Without a [V£)] mM -Bound 
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the Descent Velocity Bound 

Figure (4.6.41) compares the time variations of the thrust coefficient obtained from 
optimal control programs with and without a descent velocity bound. The latter 
consists of an initial drop in the collective pitch control. This is followed by a 
gradual, steady increase in the thrust coefficient till the final touchdown. The case 
with the descent velocity bound is characterized by phases: 

(1) The entry phase consists of a similar sharp drop in the thrust coefficient 
to preserve the rotor RPM. This is followed by a steady increase in the 
thrust until the steady state value of 0.08 is reached. 

(2) In the steady descent phase, the thrust coefficient is maintained at a 
constant value of 0.08. 

(3) The landing maneuver ends with a rapid increase in the collective pitch 
when touchdown is imminent. This rapid control movement is a special 
feature in the landing flare phase of the optimal, autorotative descent of a 
helicopter with a descent velocity bound. 

As was found in the case without a descent velocity bound, the horizontal component 
of the thrust coefficient Cj,, shown in Figure (4.6.39) is initially positive for this 
particular case with a low entry speed. The helicopter is initially accelerated to a 
higher forward velocity. The thrust vector must eventually be rotated backward to 
decelerate the vehicle for a safe landing. This reversal in the sign of the horizontal 
thrust component has been delayed when compared with results found in the case 
without the descent velocity bound. 
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The time histories of the descent velocity, forward speed, rotor RPM, as well as 
the optimal trajectory obtained for the nominal case are given in Figures (4.6.42)- 

(4.6.43) . 

The time histories of the helicopter's descent velocity obtained from optimal control 
programs with and without an upper bound on the descent velocity is compared 
in Figure (4.6.44). Those for the helicopter's forward speed and the angular speed 
of the main rotor are compared in Figure (4.6.45). It is apparent from Figure 

(4.6.44) , that “peaking’' of the descent velocity has been effectively suppressed 
after the introduction of an upper bound on the descent velocity. To cover the 
same vertical distance of 423 feet, the flight time in the case with a [Vx>] mai -bound 
has been lengthened. The touchdown time has increased from 11.3 seconds in the 
case without the !^ x>] r7jax -bound to 16.3 seconds in the case with the bound. 

The time variations of the helicopter’s forward speed with and without the descent 
velocity bound are compared in Figure (4.6.45). These variations are qualitatively 
similar. In both cases, the forward speed of the helicopter first increases, reaches a 
maximum before it is decreased to near zero at touchdown. The time at which the 
maximum forward speed is reached in the case with the descent velocity bound has 
been delayed when compares with that without the descent velocity bound. This 
is due mainly to the late rotation of the thrust vector to the rearward direction as 
explained before. 

Differences in the time variations of the main rotor RPM in cases with and without 
the lVx>j max -bound are most significant. In the case without the [V£>] m0X -bound, the 
peak angular speed of the rotor is about 30 percent higher than its nominal value 
of 353 rpm. Figure (4.6.45) clearly shows that this “peaking” in the angular speed 
of the rotor has been suppressed as a result of bounding the vertical sink-rate. The 
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Figure 4.6.42 Time Variation of the Descent Velocity and the 
Optimal Flight Profile 
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Figure 4.6.43 Time Variations of the Forward Speed and the 
Rotor RPM 
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Figure 4.6.44 Comparison of the Time Variations of the Descent Velocity 
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Figure 4.6.45 Comparison of the Time Variations of the Forward Speed 
and Rotor RPM with and without a (Vp] max -Bound 
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peak value of the angular speed in the nominal case is actually less than the nominal 
value. This suggests that the peaking of the rotor RPM found in the solution of the 
problem without a descent velocity bound can be eliminated in two different ways. 
The overspeeding of the rotor can be removed either directly by the addition of an 
inequality constraint on the rotor RPM or indirectly through the use of an upper 
bound on the descent velocity. 

4.6.3 Some Generalizations 

If we impose path inequality constraints on both the descent velocity of the heli- 
copter and the angular speed of the main rotor, the time variations of the descent 
velocity and the rotor angular speed in the optimal results might look like those 
shown in Figure (4.6.46). A plot of the angular speed of the rotor with respect to 
the descent velocity of the helicopter in the resultant optimal solution is as shown 
in Figure (4.6.46). Depending on values of [Vp] mai and flmaz, the Vp-fl plot might 
consist of one or more “corner” points. In Figure (4.6.46), the entry phase starts 
from the initial entry condition (point 0) to the point when the upper bound on 
the descent velocity is reached (point 1). Thereafter, the descent velocity is main- 
tained at a constant value while the angular speed of the rotor is increased till it 
reaches its upper bound n m0 z, at point 2. The angular speed of the rotor remains 
unchange from point 2 to 3 while the descent velocity of the helicopter is dropped. 
During this landing flare phase (from point 3 to 4), the descent velocity is being 
continuously reduced to its near zero value at touchdown. During the same period 
of time, the rotor RPM is also being reduced, indicative of the extraction of the 
stored rotational energy to cushion the landing. 

4.6.4 Effects of Perturbed Initial or Terminal Conditions 

The nominal case considered in the above subsections is with a 7.7-knots forward 
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airspeed and at an initial altitude of 423 feet. To investigate effects that initial 
and terminal conditions have on the control scheme, we consider the following three 
cases: 

(1) The initial airspeed is increased from 7.7 knots to 15 knots. 

(2) An additional terminal, horizontal distance constraint of 623 feet is added 
to the nominal problem (which is without a terminal distance constraint). 

(3) The entry altitude is increased from 423 feet to 460 feet. 

Results obtained for the first case are given in Figures (4.6.47)-(4.6.49). Those 
obtained for the second case are presented in Figures (4.6.50)-(4.6.53). Results 
obtained for the case with a change in the entry height are shown in Figures (4.6. 54)- 
(4.6.57). 

The optimal results obtained for these three cases agree qualitatively with those 
found in the nominal case. In all the cases considered, the path inequality con- 
straint on the vertical sink-rate of the helicopter has been effectively enforced by 
the optimization algorithm. Since the perturbations in initial and terminal condi- 
tions are relatively small, only minor control adjustments are needed to accomodate 
these changes. It is also interesting to note that most of these control adjustments 
are made in the entry and/or landing flare phases of the descent maneuver. Both 
the horizontal and vertical components of the thrust coefficient in the steady descent 
phase of the maneuver remain practically unchanged with the perturbed boundary 
conditions. 

The nominal problem (without a terminal, horizontal distance constraint) has a 
touchdown distance of about 560 feet. This horizontal distance has been extended 
to 635 feet in the second case. This extension in the horizontal distance has been 
achieved by an increase in the forward speed of the helicopter (see Figure (4.6.52)). 
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Initial Forward Speed = 26 fps 
Initial Forward Speed = 13 fps 


Figure 4.6.47 Comparisons of Optimal Results Obtained at Two 
Different Entry Airspeeds 
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1 Initial Forward Speed * 26 fps 

Initial Forward Speed * 13 fps 

Figure 4.6.48 Comparisons of Optimal Results Obtained at Two 
Different Entry Airspeeds 
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Initial Forward Speed = 26 fpa 

— Initial Forward Speed « 13 fps 

Figure 4.6.40 Comparisons of Optimal Results Obtained at Two 
Different Entry Airspeeds 
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— With Horizontal Distance Constraint .at 635 faat 

— • Without Horizontal Distance Constraint 

Figure 4.6.50 Comparisons of Optimal Results Obtained With and 
Without a Horizontal Distance Constraint 
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With Horizontal Distance Constraint at 635 feet 
Without Horizontal Distance Constraint 


Figure 4.0.51 Comparisons of Optimal Results Obtained With and 
Without a Horizontal Distance Constraint 
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1 1 With Horizontal Distance Constraint at 635 feet 
Without Horizontal Distance Constraint 

Figure 4.0.52 Comparisons of Optimal Results Obtained With and 
Without a Horizontal Distance Constraint 




Figure 4.0.53 Optimal Results obtained with 
a Horizontal Distance Constraint 
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Entry Height * 435 feet 
Entry Height * 460 feet 


Figure 4.6.54 Comparisons of Optimal Results obtained at Two 
Different Entry Altitudes 
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Entry Height * 435 feet 
Entry Height = 460 feet 


Figure 4.0.55 Comparisons of Optimal Results obtained at Two 
Different Entry Altitudes 
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Entry Height * 435 feet 
Entry Height * 460 feet 


Figure 4.6.56 Comparisons of Optimal Results obtained at Two 
Different Entry Altitudes 
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Figure 4.6.57 Optimal Results obtained at the 
Perturbed Entry Altitude 
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and by an increase in the overall flight time of the descent. The resultant flight 
profile is shown in Figure (4.6.53). 

A summary of the optimal results obtained in these three cases are given in Table 
(4.6.5) together with that found in the nominal case. 

4.6.5 General Conclusions 

One of the distinct features in the optimal, autorotative control of a helicopter 
with a descent velocity bound is a clear division of the landing maneuver into 
phases. These control phases, resemble those practiced by helicopter pilots in their 
autorotational training. Therefore, this control scheme should be more acceptable 
to helicopter pilots. 

A second feature of the control scheme is the suppression of the overspeed of the 
angular speed of the main rotor during the rearward cyclic flare. The increased 
flight time of the autorotational maneuver is an added benefit of putting a bound 
on the rate of descent. 

One of the major disadvantages of the control scheme considered is the high rate of 
collective pitch input observed during the landing flare phase. Another disadvantage 
of the control scheme is in the use of a large amount of cyclic pitch during the last 
few seconds of the landing maneuver. This will cause a loss of all the ground 
references during the final touchdown. However, this will not create too much of a 
problem in a practical situation since the helicopter must always be rotated to its 
level attitude before a touchdown can be made. 




INITIAL CONDITIONS! TERMINAL CONDITIONS 



Table 4.6.5 Summary of Results Obtained With and Without 
a Descent Velocity Bound 







Chapter 5 

Conclusions and 
Recommendations for Further Research 


§5.1 Summary 

5.1.1 Adequacy of a Point Maes Model in the Optimal Helicopter Landing 
Study 

A point-mass model of an OH-56A helicopter was used in the optimal helicopter 
landing study. The states were vertical and horizontal velocities, vertical and hori- 
zontal displacements, and the rotor angular speed. The cost function was a weighted 
sum of the squared horizontal and vertical components of the helicopter velocity 
at touchdown. The controls were horizontal and vertical components of the thrust 
coefficient. 

Optimal trajectories were calculated for entry conditions well within the H- V restric- 
tion curve, with the helicopter initially in hover or in forward flight. The optimal 
solutions exhibited control techniques similar to those used by helicopter pilots in 
actual autorotationai landings. The results confirm the need to drop collective pitch 
immediately after engine failure. During the landing flare phase, the thrust vector 
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is rotated to the rear in order to reduce th< forward velocity and. just before touch- 
down. the stored rotational energy in the rotor is traded for additional thrust to 
reduce the vertical velocity. The correlation between flight data and the optimal 
results establishes the adequacy of the use of a point mass model in the optimal 
helicopter landing study. 

5.1.2 The Need for Path Inequality Constraints on Both the Control and 
the State Vectors 

In order to minimize the cost function defined above, a larger rate of descent and 
a larger rotor RPM at the point of flare are desirable. However, a high rate of 
descent over a substantial period of time is unacceptable to helicopter pilots, while 
a large rotor RPM threatens the structural integrity of the rotor system. A unique 
feature of the present formulation is the addition of path inequality constraints 
on components of both the control and the state vectors. The control variable 
inequality constraint is a reflection of the limited amount of thrust that is available 
to the pilot in the autorotational maneuver without stalling the rotor. The state 
variable inequality constraint is an upper bound on either the vertical sink rate of 
the helicopter or the rotor angular speed during the descent. 

With these bounds on the control and the state vectors, the optimal solutions 
obtained will realistically reflect the limitations of the helicopter and its pilot. The 
optimal solutions consist of subarcs which are connected at suitable corners. The 
subarcs are either unconstrained, or are on the upper bound of the thrust coefficient, 
or are on the bound on the vertical sink rate. The results exhibit division of the 
landing maneuver into entry, steady descent and landing flare/touchdown phases. 
This resembles the techniques taught to helicopter pilots, and so should be more 
acceptable to them. 
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5.1.3 Comments on the Sequential Gradient Restoration Technique 

Slack variables were employed to convert the path inequality constraints men- 
tioned above into path equality constraints. The resultant two-point boundary- 
value problem with path equality constraints was solved using the Sequential Gra- 
dient Restoration Technique. The effectiveness of the SGR technique for problems 
with of both control and state variable inequality constraints was demonstrated by 
the present study. 

In general, the amount of computation involved with the use of the SGR algorithm 
is proportional to the square of the dimension of the state vector [42]. Therefore, 
the use of by auxiliary states in problems with state variable inequality constraints 
increases the computational requirements. In this regard, a transformation tech- 
nique for optimal control problems with partially linear state variable inequality 
constraints is strongly recommended [43]. The transformation technique takes ad- 
vantage of the partial linearity of the state inequality constraint so as to yield a 
transformed problem characterized by a new state vector of minimal size. Substan- 
tial savings in computer time can be achieved with this transformation technique. 

5.1.4 Reduction in the H-V Restriction Curve and Other Potential Ap- 
plications 

Even though effects of the pilot time delay and other factors have not been taken 
into account, the present study indicates that a substantial reduction might be 
achievable in the H-V restriction zone using optimal control techniques. Results 
computed using the optimal technique thus provide a benchmark for comparisons 
with other control techniques. These optimization techniques could also be used to: 

hpln inctmrt nilots on pood mit.orot.»t.ion t.prhnintw* 

\ - / r r o x * 

(2) reduce the risk/time/effort involved in establishing the H-V restriction 
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zones by flight tests. 

(3) provide an objective comparison of the autorotation capabilities of 
different helicopter models. 

(4) assess the influence of vehicle parameters on autorotation during 
preliminary design. 

§5.2 Recommendations for Further Research 

5.2.1 Refinements of the Mathematical Model 

The accuracy of studies like this one could be improved by adding the following 
refinements to the mathematical model: 

(1) induced velocity dynamics (cf. references [15], [19], [20] and [70]): 

(2) ground effect (cf. references [15 and 7lj): 

(3) variation of the profile drag coefficient with rotor blade angle of 
attack (cf. reference [24 ); 

(4) rigid body dynamics, which would show the effects of pitch 
attitude during descent and at touchdown. 

These refinements increase the dimension of the state vector (refinements (1) and 
(4)) and also the complexity of the analysis (refinements (2) and (3)). They are 
considered to be of secondary importance (cf. Section (2.2)), and should be included 
only if improvement in the accuracy of the optimal programs is required. 

5.2.2 The Effects of the Engine Failure Mode (cf. [5]) 

The nature of the power reduction transient affects pilot technique and the response 
of the helicopter. Common causes of power loss are fuel starvation, fuel control mal- 
functions, engine deterioration, and damage from an external source. The transient 
nature of the power reduction varies greatly w r ith the cause and directly influences 
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the necessary corrective pilot actions. The pr< >ent study considered the most crit- 
ical case of a sudden complete power loss, caused as a result of an engine seizure, 
disintegration, or drive system failure. It would be of interest to study other “less 
critical” cases. 

5.2.3 The Effects of Pilot Time Delay (cf. [5]) 

The term “time delay” defines the time lag between the instant of power loss, and 
the time recovery action is initiated by the pilot. It consists of the recognition time 
and the reaction time of the pilot. Therefore, the amount of time delay involved is 
a strong function of both the engine failure cues (either audio, visual or kinesthetic 
signals) available to the pilot and the experience of the pilot involved. 

A time delay on the order of 1.5-2 seconds is typically used in the establishments of 
the H-V restriction curves. It is recommended that effects of the pilot time delay on 
the autorotational landing of a helicopter from critical flight conditions be studied. 
This can be studied through the incorporation of a time delay into the helicopter 
landing problem formulated in Section (2.9). 

5.2.4 Pilot-in-the-loop Simulation 

The importance of pilot-in-the-loop simulations for helicopter research has been 
emphasized in references [4, 10, 15 and 65] . The usefulness of simulations in autoro- 
tation research has been confirmed in reference [15]. Simulations provide important 
information about pilot workload, pilot ability to track a given optimal flightpath, 
limitations of hardware involved, and the structural integrity of the vehicle. They 
may provide other information, such as whether a particular sink rate or pitch 
rate is acceptable to helicopter pilots. Data acquired from these flight simulations 
and comments from helicopter pilots should be studied carefully and appropriate 
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modifications made in the formulation of the optimal helicopter landing problem. 
Only through iterative cycles of this kind can a realistic optimal control program 
be established and implemented. 
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Appendix A 

Verification of 
Point Mass Helicopter Model 


§A.l Introduction 

A simple point mass helicopter model is used in the study of optimal autorotative 
control of the OH-58A helicopter. It includes the effects of gravity force, main rotor 
thrust, parasite drag (on the fuselage) and profile drag (on the rotor blades). The 
effects of increased induced power loss when the helicopter operates in the vortex- 
ring state have also been included. Secondary effects, such as induced velocity lag. 
ground effects etc. have been neglected. 

Experimental data describing the steady state sink rate of the standard OH-58A he- 
licopter in autorotation are given in References [12,18]. Fig.(A.l.l) (from Reference 
[18]) shows the steady autorotational sink rate of the standard OH-58A helicopter 
at a constant rotor speed of 354 rpm.. A comparison of these flight data with those 
computed by the simplified point mass model is made to establish the validity of 
the point mass model, at least in the steady state. The result of such a comparison 
is given here. 
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Figure A. 1.1 Steady Autorotational Sink Rtae of the Standard 
OH-58A Helicopter 
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§A.2 Analysis 

Table (A. 2.1) contains a list of parameters and their values used in the point mass 
model. Other than the estimated values of 24 ft 2 for the equivalent flat plate area 
of the fuselage / e , and 0.0087 for the mean profile drag coefficient of rotor blades 
all the other values were taken from Reference 18 . These two estimated values 
■were obtained by fitting the flight data for steady sink rate with the steady sink 
rate from the point mass model. 

In the steady state, the angular speed of the rotor, as well as the forward and 
vertical speeds of the helicopter are constant. We can therefore equate the first 
three of the five equations of motion derived in Chapter (2) to zero. The resulting 
equations are : 


c 0 * 0.05( 


go - mo{CT,*j 2 + fw\vu 2 -I- u 2 ) 

Ct z <^ 2 — /u\/ w 2 -|- u 2 
Ct z m - Ct z tv 


) ■+■ k 0 fj {Cj x + Cj m ) * 


0, 

0, 

0. 


where the definitions of the constants are : 


go = 


m 0 = 



*o = 


Co = 


ko = 


20000 

opT’ 

2pnR 3 

m 

5/ e 

47ri?2 ’ 

pirR s 

10 Ir ’ 

|<7Crf(l0 3 ) , 

Kind 
\/2000 ’ 


Po = 0.05(>/2000) . 
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S/NO. 


SYSTEM PARAMETER 


VALUE USED 





f a . equivalent flat plat* araa (ft 2 ) 


g. acceleration due to gravity (ft/aee 2 ) 


P. air density (a lug/ ft 3 ) 




a, mass of helicopter 


0. rotor eolidity 


ib.vt.) 








. induced velocity correction factor 


*• & (e f> 


R. radius of main rotor (ft) 


u . nominal rotor speed (rpa) 


6 m . profile drag coefficient 
* (NACA 0012 airfoil assuaed) 


7. rotor system Loek number 




92.17 


0.002378 


93.16 

(3000.0) 


0.048 


1.13 



17.63 


3S4.0 


0.0087 



Table A. 2.1 System Parameters 
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These three nonlinear algebraic equations, with their five unknown («, w, u, Cj t . 
and Cj t ). must be solved simultaneously. By fixing the rotor speed at a constant 
value of 354 rpm., then for any given value of forw ard speed, u, the above system 
admits an unique solution for it, Cj,, and Cj x . Alternatively, the airspeed of the 
helicopter can be fixed at a constant value of 45 Knots and the variation of the 
steady state sink rate with the rotor rpm. can be studied. Computer programs 
have been written to solve these simultaneous nonlinear equations. The results are 
given in Tables (A.2.2)-(A.2.5). 

§A.3 Comparison of Computed Results with Flight Data 

The effects of parasite and profile drag on the steady-state sink rate of the OH-58A 
helicopter are given in Figs. (A. 3.2) and (A. 3.3) respectively. The profound effect of 
increased parasite drag at high speed can be seen clearly in Fig.(A.3.2). Fig. (A. 3.3) 
shows the rise in the steady state sink rate needed to balance the higher profile 
power loss as the mean profile drag coefficient is increased. 

By trial and error, it was determined that a good fit between the computed and 
measured sink rate is obtained with 0.0087 and f t = 16.0 ft 2 . A comparison 

of the computed results with the flight data is given in Fig. (A. 3.4). The scatter of 
the experimental data shown in the figure is possibly due to the variation of the 
side-slip angle of the helicopter during the tests. As can be seen, the computed sink 
rate falls within the range of the experimental data. 

Further confirmation of these results can be obtained by fixing the forward speed 
at a constant value of 45 Knots, and computing the steady sink rate at different 
rotor speeds. The computed results, plotted in Fig.(A.3.5), compare well with those 
found experimentally (see aiso Fig.(A.l.l)). 
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SPEED 

( KNOTS ) 

SINK RATE 
( FPM ) 

CT ( Z ) 

CT ( X ) 

IS . 4 

2138 

0.062 

0.000 

17.4 

2018 

0.062 

0.001 

19.3 

1914 

0.062 

0.001 

21.2 

1825 

0.062 

0.001 

23.2 

1749 

0.062 

0.001 

2 S .1 

1683 

0.062 

0.001 

27.0 

1628 

0.062 

0.001 

29.0 

1581 

0.062 

0.001 

30.9 

1541 

0.062 

0.001 

32.8 

1508 

0.062 

0.001 

34.8 

1481 

0.062 

0.001 

36.7 

1459 

0.062 

0.002 

38.6 

1442 

0.062 

0.002 

40.6 

1429 

0.062 

0.002 

42 . S 

1421 

0.062 

0.002 

44.4 

1417 

0.062 

0.002 

46.4 

1417 

0.062 

0.003 

48.3 

1421 

0.062 

0.003 

50.2 

1429 

0.062 

0.003 

52.2 

1440 

0.062 

0.003 

54.1 

1454 

0.062 

0.003 

56.0 

1472 

0.062 

0.004 

58.0 

1494 

0.062 

0.004 

59.9 

1519 

0.062 

0.004 

61.9 

1548 

0.062 

0.004 

63.8 

1581 

0.062 

0.005 

65.7 

1617 

0.062 

0.005 

67.7 

1657 

0.062 

0.005 

69 6 

1700 

0.062 

0.006 

71.5 

1748 

0.061 

0.006 

73 . S 

1800 

0.061 

0.006 



0.062 

0.062 

0.062 

0.062 

0.062 

0.062 

0.062 

0.062 

0.062 

0.062 

0.062 

0.062 

0.062 

0.062 

0.062 

0.062 

0.062 

0.062 

0.062 

0.062 

0.062 

0.062 

0.062 

0.062 

0.062 

0.062 

0.062 

0.062 


Table A .2.2 Variations of Steady State Sink Rate With Forward Speed 

at Constant Rotor Speed of 354 rpm [ f e — 16 ft 7 , 6 t = 0.0087 ] 
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SPEED 

(KNOTS) 

DE-0.007 

SINK RATE(FPM) 
DE«0.0067 

DE*0.010 

0.00 

2801 

2843 

2674 

9.86 

2760 

2803 

2634 

7.73 

2632 

2679 

2713 

11.60 

2409 

2464 

2503 

15.47 

2053 

2138 

2202 

19.34 

1808 

1914 

1993 

23.21 

1632 

1749 

1637 

27.08 

1505 

1628 

1721 

30.95 

1415 

1541 

1637 

34.81 

1352 

1481 

1579 

38.68 

1311 

1442 

1542 

42.55 

1290 

1421 

1522 

46 42 

1284 

1417 

1519 

50.29 

1295 

1429 

1531 

54.16 

1319 

1454 

1557 

56.03 

1359 

1494 

1596 

61.90 

1412 

1548 

1653 

65.76 

1480 

1617 

1722 

69 63 

1562 

1701 

1807 

73.50 

1660 

1800 

1907 


Tible A. 2.3 Variation of Steady State Sink Rate With Forward Speed at 
Constant Rotor Speed of 354 rpm ( constant f t = i6 fi 2 j 



SPEED 

(KNOTS) 

FA*16.0 

SINK RATE (FPM) 
FA«24.0 

FA»32 . 

0.00 

2843 

2835 

2805 

3.86 

2803 

2795 

2787 

7.73 

2679 

2671 

2663 

11.60 

2464 

2455 

2447 

15.47 

2138 

2134 

2131 

19.34 

1914 

1916 

1918 

23.21 

1749 

1757 

1766 

27.08 

1628 

1645 

1662 

30.95 

1541 

1569 

1597 

34.81 

1481 

1521 

1563 

38.68 

1442 

1498 

1556 

42.55 

1421 

1497 

1576 

46.42 

1417 

1517 

1621 

50.29 

1429 

1556 

1691 

54.16 

1454 

1616 

1786 

58.03 

1494 

1695 

1909 

61.90 

1548 

1795 

2061 

65.76 

1617 

1917 

2246 

69.63 

1701 

2063 

2468 

73.50 

1800 

2235 

2733 


Table A. 2.4 Variation of Steady State Sink Rate With Forward Speed at 
Constant Rotor Speed of 354 rpm [ constant 6 e = 0.0087 ] 
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SPEED 

(KNOTS) 

SINK PATE 
(FPM) 

CT(Z) 

CT(X) 

CT 

300 

1155 

0.086 

0.003 

0.086 

303 

1170 

0.084 

0.003 

0.084 

307 

1185 

0.082 

0.003 

0.082 

310 

1201 

0.080 

0.003 

0.080 

314 

1217 

0.079 

0.003 

0.079 

317 

1233 

0.077 

0.003 

0.077 

321 

1250 

0.075 

0.003 

0.075 

32S 

1267 

0.074 

0.003 

0.074 

328 

1264 

0.072 

0.003 

0.072 

332 

1302 

0.070 

0.003 

0.070 

335 

1320 

0.069 

0.003 

0.069 

339 

1339 

0.067 

0.003 

0.066 

342 

1356 

0.066 

0.003 

0.066 

344 

1377 

0.065 

0.003 

0.065 

349 

1397 

0.063 

0.002 

0.063 

353 

1417 

0.062 

0.002 

0.062 

355 

1436 

0.061 

0.002 

0.061 

360 

1459 

0.060 

0.002 

0.060 

363 

1460 

0.059 

0.002 

0.059 

367 

1502 

0.057 

0.002 

0.057 

370 

1524 

0.056 

0.002 

0.056 

374 

1547 

0.055 

0.002 

0.055 

378 

1570 

0.0S4 

0.002 

0.054 

381 

1593 

0.053 

0.002 

0.053 

385 

1617 

0.052 

0.002 

0.052 

388 

1642 

0.051 

0.002 

0.051 

392 

1667 

0.050 

0.002 

0.050 

395 

1692 

0.049 

0.002 

0.049 

399 

1718 

0.049 

0.002 

0.049 

402 

1744 

0.048 

0.002 

0.046 

406 

1771 

0.047 

0.002 

0.047 


Table A. 2.5 Variation of Steady State Sink Rate With Rotor Speed 
at Constant Forward Speed of 45 Knots 
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Figure A.S.2 Effects of Parasite Drag on Steady State Sink 
Rate of Standard OH-58A Helicopter 
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2-10 


A. Verification of PoiDt Mass Helicopter Model 


ORIGINAL & 

OF. POOR QUALITY 



\ 


Figure A.S.4 A Comparison of Computed Sink Rate With 
Experimental Results of Reference [18] 
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Figure A.S.S Variation of Sink Rate With Rotor Speed at 
a Constant Forward Speed of 45 Knots 
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A. Verification of Point Mass Helicopter Model 


§A.4 Conclusions 

The Flight test data were used to validate the point mass model and also to provide 
a basis for establishing values for certain parameters in the model. The confidence 
that can be placed in the results of the optimization programs is considerably en- 
hanced. 


Appendix B 


Supporting Analysis 


§B.l Introduction 

The landingof a helicopter after engine failure has been formulated as an optimal 
control problem with a path equality constraint (cf. equations (57-64) in Section 
(2.9)). The problem is repeated here for ease of reference: 

min / = Uz XJ 2 + W,x 2f 2 ) = *(*,), 

u, * * 

^ = (*1 *2 *3 *4 *l) T , 

H = (ill «2 «3) r » 
f=(r/). 

where 0, jt and 9 are the control .state and parameter vectors. 


subjected to : 
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B. Supporting Analysis 


(1) equations of motion ( X v — / )’• 

xi v = T f {go - mo(tii*3 2 + /* i\/*i 2 + *2 2 )), 
x 2 v = r/mo(u2*3 2 - /*2\/*i 2 + *2 2 )> 

Z 3 V = — T/t'oX3 2 (CO + A\/ui 2 + U 2 2 ), 
x 4 V = O.lr/Xi, 

X5 V = 0.1r/X 2 . 

The initial condition of X is given by: 

X 0 = (0,u 0 ,l,0,0) r . 

(2) path equality constraint (S{X,U t ir) =0): 


(«i 


> a =°- 

Lj, 


(3) terminal constraints {\p{Xf,x) = 0): 


Xif — hf — 0, 

*3/ - df = 0. 

Id this appendix, detailed expressions for (5x3 matrix), (5x5 matrix), §£ 
(5x1 matrix), |§ (1x3 matrix), (1x5 matrix) and (2x5 matrix) are given. 


§B.2 Determination of (//) f 

The induced velocity parameter // is a function of *1 and *2 which in turn are 
funtions of «i, v 2 , * 1 , * 2 . and * 3 - fi is thus a function of these variables and the 
partial derivatives of // with respect to u i,u 2 ,xi,x 2 and X3 are given below. 


B.2 Determination of (//){ 
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Let us first define 0 and 7 as follows: 


0 = *1 = Po 


7 = *2 = PO 


X2«2 ~ 

X$(tti 2 + U 2 2 ) 
«2«l + Xl«2 

*s(«l 2 + «*2 2 ) 



(i) 

4 


T’ 

(2) 


Once again, the value of the induced velocity parameter //, is determined from the 
expressions: 

f 1(0,1) = ( 1 *°/\/('/ 2 + (0 + //)*)» M ( 2 0 + 3) 2 + t 2 > 1; (3) 

| 0(0.3730* + 0.5987 2 - 1.991), otherwise. 


In the region where the momentum theory is applicable, (i.e. when (20 + 3) 2 4*t 2 > 
1) expressions for //^ and // 7 are given by the following expressions: 


where 


ftfi = 


l-z’'0 + fj' 

(0 ± ft) 

( 7 2 + (0 + fl?V 


= -(/? + //)//*• 


(4) 

(5) 

(6) 


In the vortex-ring state, the expressions of //^ and // 7 are given simply by: 


f J0 - I.U90 7 + 0.5987 2 - 1.991, (7) 

f h = l.\960r (8) 


We can then use the chain-rule: 


(//)( ) = flfi(0)( ) + ) 


(9) 
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where ( ) stands for uj, «2» *i» *2 or **• To determine (//)( j, expressions of 
j and 7( ) are therefore needed. These expressions, which are identical for both 
the vortex-ring and momentum theory states, are as given below: 


/?*, 



( 22!±1 - ) p 

X3(«i 2 + U2 2 )« 

+ f E^ 2 ) 

X3(ui 2 + U2 2 )« 


( 


( 


P o 


X 3 (U ! 2 + U 2 2 )< 

P 0 

X3 ( u l 2 + «2 2 )« 


)(— xi - 1.5«i 
)(+X2 — 1.5«2 


(x 2 tt 2 -g|Ul) 

(til 2 + M2 2 ) 
(x 2 tt2-XiUi) 
(tti 2 + U2 2 ) 


). 


), 


similarly the expressions for ^ j are given by: 


( 10 ) 


7*i 

7*, 

7*j 

7.i 

7., 


+{ ) 

X3(uj 2 + tt2 2 )< 

+( ) 

X3(ui 2 + U2 2 )< 


( ^ 

X 3 (tij 2 + U2 2 )< 

( £2 

X3(Ui 2 + «2 2 )< 


)(+X 2 - 1.5«i 
)(+Xi - 1.5«2 


(x 2 m +X!U 2 ) 

(Ul 2 + «2 2 ) 

(x 2 tti -4- x t u 2 ) 

(tti 2 + U 2 2 ) 


). 


), 


(H) 


Expressions of (//),, , (//) XJ , (//) XJ , (//)., and (//). a can therefore be computed 
using equations (4-11) for either the vortex-ring or momentum theory state. These 
results will be needed in equations (12) and (13) in the next two sections. 


§B.3 Expressions for df/d& , dfjdit and df/d x 


B.3 Expressions for df/dU_ , djjdi and Of/dif 
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^ is a (5x3) , a (5x5) and §£ a (5x1) matrix. Detailed expressions of these 
matrices are as given in the following sub-sections. An individual element of a 
matrix is indicated (in the usual way) by its row and column position in the matrix. 
An element which is identically zero is not given. 


B.3.1 Expressions for df/dO . 


is a (5x3) matrix. Non-zero elements of are given below: 


(1.1) = -T f m 0 xz 2 t 

(2.2) = +T f moxz 2 , 

(3, 1) = — T/ioX 3 2 (- ~~ + 1.5Jko//( 


( 12 ) 




Xi («l 2 + «2 2 )^ 

+ 1.5ho//(— ■ — ) + ko{ «i 2 + u 2 2 ) *(/»„,). 

(u 1 2 + u 2 2 )« 


(3,2) = - Tf i 0 x 3 2 (+^± - ’ " “ 2 

*3 


B.3. 2 Expressions for df/dX . 

f£ is a (5x5) matrix. Non-zero elements of are given below: 


(1, 1 ) = -m 0 r//(\/r i 2 + x 2 2 + — p== = ), 

V*i 2 + x 2 2 


(1,2) = -m 0 r//(- 


*1*2 


), 


V*l 2 + *2 2 

(1,3) = — 2T/motiiX3, 


(2.2) = -m 0 f / /(N/r|i + ^+ *1, , ), 

V*l 2 + *2 2 

(2.3) = 2r^mott2X3, 

(3. 1) = -r / * 0 * 3 2 (_2^Iii + ko(fi) Ml ( ui 2 + u 2 2 )* ), 

*3 

4ft 0.01 u* « . ) 

(3.2) = -r/ioxj'l+^i + *»(//)., (“i* + OJ 1 )*), 

*3 


( 13 ) 
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(3,3) = -T/io(2co* 3 + 0.01 (« 2 «2 - *i«l) + *0*3(«i 2 + «2 2 )« (2// +*3(//),,)), 

(4.1) = 0.1f/, 

(5.2) =0.1r/. 


B.3.3 Expressions for df /dr. 

is a (5x1 matrix). Non-zero elements of §£ are given below: 


(1. 1) = go - m 0 (ttiz 3 2 + /xiV^*i 2 + * 2 2 ), 

(2. 1) = m 0 (u 2 i 3 2 - /* 2 \/*i 2 + *2 2 ), 

(3. 1) = -i 0 z 3 2 (c 0 + A\Zui 2 + u 2 2 ), (15) 

(4.1) =0.1zi, 

(5.1) =0.1*2- 


Expressions for the inflow ratio (A) and advance ratio ( fi) are given by the following 


expressions: 


A = 0.01 ( 

/i = 0 . 01 ( 


X 2 U 2 - *lU t 

*3\/ui 2 + U 2 2 
*2Uj -f *ltt2 

*3\/t*l 2 + «2 2 


) + ^0//(«l 2 + «2 2 ) 
)• 


1 
« . 


( 16 ) 


The first equation of (16) on A is needed in equation (15). The second expression 
on fi will be needed in equation (19) of Section (B.6). 


§B.4 Expressions for dS/dU 

The path equality constraint 5 is not a function of either £ or f but only a function 
of U. Expressions for and ff are therefore identically zero. is a (1x3) matrix. 
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Non-zero element of are as given below: 


( 1 , 1 ) = 

(1.2) = 2« 2 (i - 

(1.3) = -2« 3 . 


^T. 


(17) 


§B.5 Expressions for d<f>/dXf and drf/dj^f 

The terminal cost ( <f > ) (a scalar) and the hard terminal constraints (if>) (qxl vector) 
are usually combined to form a ((q+l)xl) vector $. $ is not a function of either 
Uotx but only a function of X f . || and |f are therefore identically zero, is 
a (3x5) matrix. Non-zero elements of are given below: 

(1.1) = *1, 

(1.2) = W t x 2 , 

(18) 

(2.4) = 1, 

(3.5) = 1. 


§B .6 Determination of Collective pitch setting, £75 

The collective pitch setting required to obtain a given amount of thrust may be 
obtained from blade element theory as (24]: 

_ (1 + §,.»)(*£) + |A (1 - 
<1 -„»+!„*) 

where £75 is the rotor collective pitch angle at 75 percent chord, while 0 and a are 
the rotor solidity ratio and rotor blade two dimensional lift curve slope respectively. 



250 


B. Supporting Analysis 


The quantities, n and A were defined by equations (16) of section (B.3.3). The 
collective pitch angle 0 7 s is therefore given by the following expression: 



(l + V) («jL!0^!iiZ) + | (l -4) A 

l- It 2 


( 19 ) 


Appendix C 

Example Problems 
Solved Using the CPF Algorithm 


In the following section, three examples are described, all pertaining to problems 
with an unknown parameter. Problem (l) is a minimum time and energy problem 
with an unspecified terminal time. Problem (2) is a specified control law prob- 
lem and Problem (3) is the classical Brachistochrone problem. These problems 
are solved numerically using the Combined Parameter and Function Optimization 
Algorithm developed in Section (3.2). 

§C.l Example Problems 

PROBLEM (1) MINIMUM TIME/ENERGY CONTROL 

Consider the following minimum time/energy problem: 

1 f tf 9 

min/ = tj° + -0 I u 2 dt 
1 2 Jo 

The equation of motion is given by 


x-u 
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C. Example Problems Solved Using the CPF Algorithm 


The initial condition is 


x(0) = 1. 

The terminal condition is 


x(t f ) = 0. 

Values of a and 3 are chosen to be 

Oc = 1 , 

3=1. 

Initial guesses for u(/) and tj are 

u(t ) = 0, for 0 < t < tf 

<f= 1- 

The problem is first converted to one with fixed end-time (cf. equation (38) of 
Section (3.2)). The transformed problem with a=0= 1 is 

If 1 2j 

mm I = tt + -tf I u ar 

1 2 1 J 0 

The independent variable is now r which varies from 0 to 1. 

The modified equation of motion is given by 

X = tfU , 

where ( ) now denotes time differentiation with respect to r. 

The initial condition of x(0) remains unchanged. 


C.l Example Problems 
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The transformed problem is one with an unknown parameter t j. The problem can 
be solved using the Combined Parameter and Function optimization algorithm and 
the optimum solution found is 

u{t) = -25 , 
tf = 2"5 . 

The minimum value of the cost function is given by 

Imtn — 25 = 1.4141 . 

PROBLEM (2) SPECIFIED CONTROL LAW PROBLEM 

Consider the following optimization problem: 

(xj 2 + u 2 )dt 

subject to the following equations of motion 

X\ = *2 + 0.01X2 3 , 

X 2 = — 4xj — 5x2 + 4u . 

Initial conditions are 

xi (0) = 2 , 

x 2 (0) = 0. 

There is no terminal condition. 

The control law is specified to have the following form 

u = 6xj , 


nin I=\ I 
u 2 J 0 


where b is an unknown constant whose value is to be optimally selected. The 
transformed problem becomes 
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C. Example Problem? Solved Using the CrF Algorithm 


If 1 

min I = - I *i ,l ■+ b 2 \di 

b 2 J 0 

with the following equations of motion 

*1 = *2 + O.OlX® , 

*2 = 4xi(5 - 1) - 5x2 . 

while the initial conditions of xi(0) and X 2 ( 0 ) remain unchanged. 

The transformed problem is a problem with an unknown parameter b but with- 
out any control. The problem can be solved using the Combined Parameter and 
Function optimization algorithm and the optimum value of b found is: 

b = -0.23387 


The minimum valut of the cost function determined is: 

hoc = 1.14284. 


Note that the value of I soc is (slightly) larger than the minimum value of the original 
problem, w’here u(t) is open, which is: 


1.13289. 


PROBLEM (3) BRACHISTOCHRONE PROBLEM 

Consider the Classical Brachistchrone Problem 


min I = tf 
9 1 


The equations of motion are given by 


C.l Example Problems 
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X] =■ X 3 COS0, 

±2 = X3 sin 6 . 
X3 = sin 6 . 

The initial conditions are given by 

xi (0) = 0, 
x 2 (0) = 0, 
* 3 ( 0 ) = 0 . 

The terminal condition is given by 


xi(</) = 1. 

This problem with an unspecified terminal time can be transformed into one with 
t f as the unknown parameter (cf. Problem (1)). The transformed problem can be 
solved using the Combined Parameter and Function Optimization Algorithm. The 
minimum value of the cost function determined is 

4mn = 1-7724 


which is very close to the exact solution of v / 7r. 



Appendix D 

Generalized Transversality Condition 


§D.l Introduction 


In the classical calculus of variations problem, we seek to minimize a cost function 
J subject to a set of differential equations and terminal constraints: problem: 

f*s 

min J = <t>{xj, tj) — / L(x.u,1)d1, 

t Jo 

x = /(x, 

x(0) = given, 

and t = 0. 


Here x(n x 1), ii(m x 1), and n[p x 1) are the state, control and parameter vectors. 
%j){q x 1) are the terminal constraint functions, tj is the specified or unspecified 
end-time of the problem. If if is unspecified, first-order necessary conditions of 
the problem are given by the following relations: 


D.l Introduction 
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X = f{x,u,t ) , 

*(0) = given , 

X T = -H r , 

H u = 0, 

xl’iifjf) = 0. 

and (-37 ■” #) = 0. 

ot t=t f 

where for convenience, we have defined scalar functions H (the Hamiltonian func- 
tion) and $ as : 


H(x.H,Lt) = L(x. u,t) + X T f(x,u,t), 

$(Xf,tf,fI) = + fJ. T ll>(Xf,tf) . 

The last condition in (1) is the so called “ transversality condition ” in the classical 
literature. It determines the value of the unspecified terminal time, tf. 

The above problem with unspecified end-time can be converted into one with fixed 
end-time with the following change of independent variable: 

t 

*/ ’ 

The transformed problem becomes: 

min J = <f>i -+- f tfL(x,u, tfT)dr, 

J 0 

(f)' = tfT ), 

x{0) = given, 

0(*/,*/) ~ °- 


and 
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D. Generalized Transversality Condition 


Here ( )' denotes differentiation with respect to the independent variable r. 

The characteristics of the transformed problem are that the end-time is now fixed 
at r = 1 but we have introduced into the problem an unknown parameter t j that 
is to be ‘"optimally” selected. 

The first-order necessary conditions of the transformed problem are similar to those 
given for the original problem. If we define the following “transformed” functions: 

H(x, u,\',tf,r) = tfL(x, u, t f r) + {\') T tff(x,u, t f r), 

= t f [I(x, U, t fT ) + (A') r /(f, U, tfT )] , 

= <p{xf,tf) + {p) T xs(x f J f ). 

Here A'(r) and /T' are the Lagrange multiplier function and Lagrange multiplier 
associated with the transformed problem. 

The necessary conditions are given by the following relations: 

(*)' = *//(*»“’ */ r )’ 
x(0) = given , 

(A'i) T = (S>x)i , (2) 

H u = 0 , 

= o, 

and (4^)1 + f Htjdr = 0. 


In this appendix we prove that necessary conditions (1) and (2) are equivalent to 
each other and that the last equation in (2) is a generalized transversality condition. 


D.2 Proof 
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§D.2 Proof 

The proof consists of four parts. The first two parts of the proof establish the 
facts that A' = A and jl' = fl. The last two parts prove that in the special case 
when x = if. the last equation of necessary conditions (2) becomes the classical 
transversality condition. 

Part (1) 

From the optimality condition of (2), we have: 

Q 

— H(x,u,X',t f ,r) 

d - T 

tf — \L{x,U,Ttf) + (X') f{x,U,Ttf )] 



Since t f is non-zero, we must have: 

d , 




du 


L{x, u, Tt f) -r (A') f(x,u,Tt f ) = 0 


( 3 ) 


From the optimality condition of the original problem, we have the following rela- 
tion: 


H u = 0 , 

— [I(x,«,t) + A 7 7(z,u,t)] = 0. 

A comparison of equation (3) and (4) leads to the conclusion that: 

A' = A 

The same conclusion can be reached by comparing the adjoint equation: 

A r = - H z , 

A 

with 


( 4 ) 


(4a) 




Part (2) 
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D. Generalized Transversality Condition 


To prove that we can compare the terminal conditions of the adjoint 

equations (i.e. the A-equation and the A'-equation). From Part (1) of the proof, we 
have: 


(A'), = (A,)' , 

- ( ar )! ’ 




( 5 ) 


Whereas the terminal condition of the A-equation is given by 


* T (1 = t f ) = \-?-(6(x f Jf) ~ fj T v{ff.1 f ))] , 
ax t=tf 


( 6 ) 


a comparison of equation (5) and (6) leads to the conclusion that 


V = V 


(6a) 


Part (3) 

Here, we prove an intermediate result that is needed in part (4). To express /J Ht,dr 
in terms of H, we note that: 

r 1 . f l d 

I Htjdr — I —\t f H(x,u,Tt f )]dT, 

f i d 

= / \H(x,u,Ttj) ■+• tf — H(x,u,rt f ) dr , 

Jo 

= J Q [ H{x,u,rt f ) + t f ^{H{x,u,t))-^-]dr. 


But since: 


t 

T ~ * 7 ’ 

dt _ 
dtf ~ T ' 

dr = — . 

l f 


and 


D.3 Conclusion 
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Therefore the above equation can be simplified to: 

■1 rt, 






if Jo 


1 f** 8 

m T, L a,W- 


u 

= j-my . 

= (H)tztj . 


( 7 ) 


This result is needed in the next part of the proof. 

Part (4) 


To express in terms of (fjf )t=t f . we note that: 

di> 


= IftjM*/'*/) “ ^ V(x/.//)).r = l , 

- t-gr )«-*/• 


( 8 ) 


From equations (7) and (8), we obtain the final result: 

(^iy)i + f (Htj)dr , 

Jo 

= l 

8t 


d9 

= <S + *)«-*, 7 


= 0 . 


(®) 


which completes the proof. 
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D. Generalized Transversality Condition 


§D.3 Conclusion 


In section (D.2). we have proved that the necessary conditions (2) for the genera) 
problem with fixed end-time (r=l) but with an unknown parameter vector (f) 
are equivalent to the necessary conditions (1) for problems without the unknown 
parameter vector tt but with an unspecified end-time //. The necessary condition 


in equation (2) that determines the value of rr is in effect a more general form of 
the classical “ transversality condition ” 



Appendix E 

Definitions 
of Matrices used in Section (3.4) 


§E.l Introduction 

In Section (3.4), we formulated a linear, time-varying Two-Point Boundary- Value 
Problem with integral path constraints. The problem was formulated in terms of 
the matrices 5, (where t = 1,...,13), A tJ and C tJ where i = 1,2 and j = 1, 2 and 
3. These matrices themselves are functions of H UX1 . H u ». H zz , H uz . H rz and H ru . 
Detailed expressions of 5,, A tJ , B tJ as well as those of H uu , H ur , H zz . H uz , H TZ 
and H nu are given in this appendix. 


§E.2 Expressions for ) 


Expressions for H UU i Hu*, H xx , H uz , H , rx and H xu are given in this section. For 
brevity, we assume that the integral payoff L is zero. Cases when L is not zero 
can be treated in either of the following ways. The problem may be converted to 
one with only terminal payoff by the addition of an auxiliary state whose integral 

r au si r / ^ T ^ 7 /_ w r / ^ 

ia aj . Aiicuiauvci), uic uiau itca j^uu\ 9tl A rtl Ji A F ) **ZX\ 9i A 9i Ji ^U7V"* A 91 j • 

L* *(p x n) and L xu (p x m) must be computed and added to the the appropriate j 
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E. Definitions of Matrices used in Section (3.4) 


matrix. We have also assumed here that the path constraints S are not functions 
of 7 t. General expression of H[ j is given by the following expression [29j: 



Without the integral payoff L , the Hamiltonian function H is given by the following 
expression: 

H = \ T f - p T S 



ORIGINAL PAGE IS 
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E.3 Expressions for S( j 


Hxiz — foifu A ">■ $U p) 


( A, dui /*, + 23 P) du^dn 


fj Jl. v>. 9 *s 3 

\ 2s I dumdXi 2s pj dumdzi 

V|=l 7=1 


21 A ia«jdi„ + 2 pJdt\Pz u ^ 

>=1 


1=1 


V \ &~fx , d S, 

1 -* * dUmOZrx ■ — * Pj dumdln 

*=1 J = 1 


Hxzip x n) 


H„ = 


d t 

gj(/jA), 

/ Ea. 5 ^ 


1=1 


1 dZjdTTj 


. is td Z} dn p 
V i=l 


2 -* A *dz n dn 1 


2/ \ 


*=1 




i=l 


if nu (p x m) 


H zz (n x n) 


#*«, = (IT.,) 3 


■ff*x - ^(/f A + 5^p) , 


symmetry 


i>.0 + £>,0 - 

*=1 t=l 1=1 t=l 


i=i 


n -o - r « 

— a - ^ , 


£A,£* + X>,§r n - , 

i=l 7=1 / 
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E. Definitions of Matrices used in Section (3.4) 


Here the time-varying functions A(r), p(r), ^{f u ) T . ^(Su) 7 ", f u ) T • jj(fu) T , 

£(S U ) 7 ’, and £(f , r ) T are evaluated on the nominal optimal path. 

§E.S Expressions for S( j 

The dimensions and detailed expressions of 5,, where t = 1 ,...,13 (note that S$ is 
not defined) are given by the following relations: 


S l (r x n) = 
S 2 {r xp) = 
S 3 (r x n) = 
S 4 (r x r) = 
5 5 (p x p) = 
Se(m x n) = 
S; (m x n) = 
5 e (m x p) = 
Sio(p x 1) = 
5n(p x p) = 
5 J2 (p x n) = 
and Sj 3 (p x n) — 


S U H 


-l 

UU 




5,, 


SuH^Hu, - 5 , . 

C fT* 

y u * 


C IJ-l cT* 

( ^ 7T 7T ) 1 + / Hjcitdr , 

Jo 

-HZUfl -Si S 4 - 1 S 3 ), 
-^ tt - u 1 (/f tt ,-5j5 4 - 1 S 2 ), 

(*ir*)l(6z)l + (t^)l(<5/7)l , 

S 5 + [ [H* u Ss - S^S 4 - 1 S 2 ]dr, 
Jo 

H„ + H wu S 6 -STs 4 - l Si, 


H ru s 7 + fi - sZs 4 -'s 3 . 


Note that 5i to S 4 are defined only when the matrix J/ uu is non-singular. This 
is referred to as the convexity condition (or strengthened Legendre- Clebsch 
condition) in calculus of variations. The existence of the matrices S$ to 5s and 
5n to S \3 also depends on the non-singularity of both H uu and S 4 (= S u HZu^I )• 
The latter condition is in effect a modified form of the convexity condition for 
problems with path equality constraints 5. It requires that all the path equality 
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E.5 Expressions for C t] 

constraints S must involve some components of the control vector u so as to ensure 
that det(S 4 ) 0. 

§E.4 Expressions for A tJ 

In terms of the matrices 5, defined in the last section, expressions for A X} are given 
below : 


A\\ (m x n) = S 6 , 

A n (m x n) = S 7 , 

An{r x n) = -S 4 ~ l Si , 

M\{r x n) = -S 4 -1 5 3 , 

Ml (m x p) = Sg , 
and A 2 z{r xp)= -S 4 ~ l S 2 - 

These expressions for A x} are needed in the computations of the neighboring optimal 
feedback gains, A, (where 1 = 1, 2 and 3) in Section (3.4). 

§E.5 Expressions for C tJ 

In terms of the matrices 5, defined in section (E.3), expressions of C t j are given by 
the following relations: 

Cn(n x n) = f x + / U S 6 , 

Cn(n x n) = f u S 7 , 

Ci 3 (n x p) = /* + f u Sg , 

Czi{n x n) = -H zx — H XU S e S^S 4 , 

C 22 (n x n) = -f x - H xu S 7 + S^S 4 ~ l S 2 , 
and C 2 3 (n x p) = —H x1t - H zu Sg -f Sj S 4 ~* S 2 . 
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E. Definitions of Matrices used in Section (3.4) 


These expression* for C t] are needed in the solution of time- varying Two-Point 
Boundary- Value Problem formulated in Section (3.4). 


Appendix F 

Iterative Solution for TPBVPs 
with Integral Path Constraints 


In Section (3.4), we consider small deviations from a nominal path produced by 
small perturbations in 6 z(0) and 6 ip. These perturbations give rise to 6x(r), 5A(r), 
$u(r), 6p[r) and 6$. The perturbed quantities are determined by a linear, time- 
varying Two-Point Boundary- Value Problem (TPBVP) with integral path equality 
constraints. The general form of the problem is given by equations (3-6) of Section 
(3.4). 

When S and /* are not functions of f , the TPBVP can be simplified and solved 

using the Backward Sweep Method of reference [29]. The resultant neighboring 
optimal feedback law has the form 


6u = -Ai 62 - \26if — A3 6*. 


The values of 6 if needed in the above feedback law can be obtained from the solution 
of integral path equality constraints. An iterative procedure which imbeds the 
solution of 6 if in the solution of the TPBVP is described in this appendix. 
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F. Iterative Solution for TPBVPs with Integral Path Constraints 


§F.l Problem Formulation 

The TPBVP consists of the following differential relations: 

© - (a %< a) (i!) <" 

and path equality constraints 

0 = 5n£* + f = 0 • (2) 

J o 

The end-conditions of (1) are 

$x(0) = given , 

(ttf)i ~ = given, (3) 

(5A)i s ($,,)i(5x)i + (jJUSji. 

Here £x(r)(n x 1), 6\{r)(n x 1) and 6n{p x 1) are variations in the state x(r), 
Lagrange multiplier A(t), and the unknown parameter ir from their nominal values. 
The matrices C, ; (r) and 5,(r) are defined in Appendix (E). 

Solution of the TPBVP with path equality constraints can be obtained using a 
modified form of the backward sweep method. The solution of the problem is given 
by the following relations: 


6X = T 6x + R6ji + P6& » 
6$ = R T 6z + Q6ji+ H6*. 


( 4 ) 


where T(n x n), R(n xq), P(nxp), Q(qx q) and H{qxp) are matrices that satisfy 
the following differential equations 


F.l Problem Formulation 
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t = -TCn - Cfi T - TC n T + C 2 1 , 

J? = -(CE + TCm)*, 

P = C23 - Ch P — T(Ci2i* + Cis) , 
Q=-R T C l 2 R , 
and H — -R T (C\ 2 P + C13 ) . 

and with the following terminal conditions 

T(l) = (*„)i , 

*(D= Wlh. 

P(i) = (*»h = o, 

<?(!) = 0, 

and H{ 1 ) = (tM 1 = 0 . 


Using these terminal conditions, the matrix time histories can be obtained by back- 
ward integrations of the differential equations given above. The task remains to 
find the value of Sir that also satisfies the integral path equality constraints 


ft = SuSx + f (Si 2 Sx + SijiA ) dr — 0. 

Jo 

Once determined, Sir can be used in the neighboring optimum feedback law 


<5u = -Ai$*- \ 2 Sf- 


( 5 ) 
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where the time*varying gains, A, ( » = 1, 2 and 3) are given by the following relations: 

Ai(m x n) = -[An + Au{T — R Q“ l Ji r )) , 

Ai{m x q) = —AnRQ~ l t 
and A 3 (m xp) = — + Mi{P — R <J“ l iT)J . 

§F.2 Solution method 

An iterative procedure can be constructed for the solution of the problem posed in 
the last section. We first note from equation (4) that 

6 A = [T - RQ~ l R T \6x + [RQ’ l \6$ + [P - RQ~ x H\6x . 

i 

Therefore the initial conditions needed for the forward integration of equation (1) 
are given by the following relations: 

(Sx) 0 = given , 

(6) 

(«)„ = [r - ft<?“'.ft r l 0 («x| 0 + [SQ-'lo^ + IP - RQ-'H^Si. 

I 

If we denote the value of Sir at the i th iteration by 5ir,, we can substitute this 
value of Sx into equation (6) and integrate equation (1) forward to obtain time 
histories of (5i(r)], and [5A(r)] t - of the i th iteration. These time histories can then 
be substituted into equation (2) to obtain the value of fl at i th iteration 

n, = 5n[5#li+ f l (Sa(r)mr)\i + Sair)[SXit))i) dr . 

Jo 

I If all the components of n, is less than e, i.e. | [H.Jy |< e, where j = l,...,p and e 

is a small, pre-selected quantity (for example e = 10" 4 ), then the value of 6x used 
in the present iteration is close enough to the correct value and can be used in the 
feedback law. Otherwise, \6x\i must be improved to better satisfy path equality 
constraints in the next (i.e. (» + l) tA ) iteration. 


! 


* 7 / 

F.2 Solution method _ 

To construct an iterative procedure, let us define the following vectors/matrices: 

© • 

P = ($*), 

A (C IX C 12 \ 

A -\C 21 C n )' 

„ _ ( <**>• 

C ~ ^ [r - RQ-'R T M6i)<, + m-'hW ) ' 

D= ([P-SQ-'^lo) 

Equations (1) and (6) can now be written in more compact forms as follow: 


jf(r) = i4(r)y + B(r)p, 
y(0) = C + Dp. 

from which we note that p (=5ir) enters into both the differential equation and its 
initial condition linearly. 

Now consider the effect of variation in p, i.e. 


p -+ p + A p, 


(where A p is not small) on the solution of equation (8). Let the corresponding 


change in y be 


y -♦ y + Ay. 


If we substitute the last two equations into equation (8), we obtain the following 
relation between the A-quantities 


Ay = /(r)Ap, 
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where flr) is a (2» x p) matrix that satisfies the following first-order, linear, time- 
varying differential equation 


/ = A{r)f + B{t) , 


( 10 ) 


and with the following initial condition 

/( 0) = [D\. 

The result given in equation (9) is needed in our subsequent analysis. 


§F.3 Iterative Procedure 

If we denote the value of f at the i th iteration by ft, and if ft does not meet the 
path equality constraints, then 


n, = SuPi + [ {Si 2 > Su) fidr # 0 , 

Jo 


To come closer to satisfying the constraints, we change p< by Api in our next 
iteration. The corresponding change in ft is given by equation (9). The path 
! constraints at the (« + \) lk iteration are 

0=n i ^=Sn(r-+AP.) + f\s l 2,S u )Wi + nr)Ap-\ir. (*» 

i 

I 

where Aft is selected in such a way that path equality constraints at the (t + l)‘ fc 
are satisfied exactly. 


If we substract the ft +1 - equation from the ft- equation, we get the following 
relation that can be used to solve for Ap». 

-n, = [Sn + J (5i2,5i 3 )/(r)<ir](Aft) , (12 ) 

therefore A ft = -K W< . 


I 


F.3 Iterative Procedure 
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where 

K = \Su+f W Sis)/ M*]-‘ 
Jo 


(13) 


Since 5 u (r), 5 i 2 (t), and Si 3 (r) can be computed (cf. Appendix (E)) and /(r) is 
the solution of equation (10), the (p x p) gain matrix K can be computed. 


Once the value of Ap, is computed, the value of p for the 
computed as 

P.+ 1 = Pi + Ap - , 

= Pi ~ Kfti • 


(» + l) 1 * iteration can be 

(14) 


The path constraints are satisfied “exactly” if we use this value of p for the (»+ l) ,fc 
iteration. 


Alternatively, an iterative procedure can be constructed as follow. Here instead of 
finding the value of K from equation (13) and using it in the computation of p ,+ 1 
for the (» + l) tfc iteration, we choose to compute the value of pi + 1 as follow: 


(p)»+i — Pi ~ ^1^1* » 
or (5ir) l+ i = (5jr) t - — X[fl]i . 


(13) 


Equation (13) states that the correction that has to be made after the i lh iteration 
is proportional to the negated value of the constraint violation ft,-. K(p x p) is 
an arbritarily selected, positive gain matrix. Now instead of meeting the path 
constraints in one single step (with substantial amount of calculation in that one 
step), we meet them iteratively in several steps. The iterative procedure ends when 
the conditions | [0,]/ |< e for; = 1.....P are met (this also means that the correction 
to p in the next iteration is going to be insignificantly small). A flow chart of this 
iterative procedure for the case when p * 1 is shown in Figure (F.3.1). 
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F. Iterative Solution for TPBVPs with Integral Path Constraints 


Once determined, the value of can be used in the neighboring optimal feedback 
law and the TPBVP with path equality constraints is solved. 


F.3 Iterative Procedure 
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Output the value 
of (5 t), 


I 


Obtain feedback lav* 
from eq. ( 5 ) 


^Sto 


PigtlTt 7f eSel 


An iterative procedure to find Sx 







Appendix G 


Example Problems 
Solved Using the SGR Technique 


Four example applications of the Sequential Gradient Restoration Optimization 
Algorithm are given in this chapter. Example (1) has a 6tate variable inequality 
constraint. Example (2) has a bound on the control. Example (3) is one with 
multiple path constraints involving both the state and the control vectors. Lasth, 
Example (4) involves a situation where there is a bound on the time rate of change 
of one of the states. Whenever necessary, the path inequality constraints are first 
converted to equality constraints that involve the controls. The Sequential Gradient 
Restoration (SGR) method is then used to solve the optimal control problem with 
nondifferential constraints. 


§G.l BOUNDED BRACHISTOCHRONE PROBLEM [29] 


Consider the following optimization problem: 


min I = 

U 
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The equations of motion are given by 

V = g sin-), 
i — V cos -y . 
y = Usin-y, 


The initial conditions are : 


U(0) = 0, 

x(0) = 0, 

y(o) = o. 

The terminal condition is : 

x{t f ) = 1. 


The minimum time problem is further subjected to a state variable inequality con- 
straint of the following form : 



The above problem is the same as the one posed in Chapter (3) of Reference [29' 
(pps 119-120) with h = 0.2 and tan# = 

The inequality constraint can be converted into a path constraint with the use of a 
Valentine-type auxiliary state z 

+ 0.2 - y - z 2 = 0. 

The initial condition of the auxiliary state z is given by 


r(0) = +\/a2 . 
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G. Example Problems Solved Using the SGR Technique 


If we take the first time derivative of the path constraint and substitute for x the 
dynamic equation of x, we obtain a path equality constraint which involves the 
control 

^ V cos *y - V sin q — 2zw = 0 , 

w’here the auxilliary control tr is defined as the time derivative of the auxiliary state 

z 

i = u< . 

Therefore the transformed problem is given by 

x 1 = t fV cos q , 

y’ = tjV sin 7 , 

z 1 = tfXV , 

V' = tfgsin~/. 

Since the end-time ofthe problem is not specified, we have made a change in the 
independent variable from t to r where 

t 

T ~ */' 

The “new" independent variable of the problem is now r (which varies from 0 to 
1 ); and ( )' denotes differentiation with respect to r. 

The initial conditions of the transformed problem are 

x(0) = 0, 

y(o) = 0, 

x(0) = +\/02, 

V(0) = 0, 

and the terminal condition is 


*(1) = I- 


G.2 BOUNDED CONTROL OF DOUBLE INTEGRATOR PLANT |29: 
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The path equality constraint of the transformed problem is 

^ V cos ~i — V sin *7 — 2zw = 0 . 

The minimum-time of this problem for the special case when 5 = 1 and / = 1 is 
found to be 1.7741. This minimum time compares well with that computed using 
the formula given in .29 (t j = 1.7795). 

The optimal trajectory of the above problem is given in Figure (G.1.1). 


§G.2 bounded control of double integrator plant 

[29] 


min I = tt 

u J 

The equations of motion are given by 

i = y, 

y = u. 

The initial conditions are : 

x(0) = 0, 

y(o) = o. 

The terminal conditions are : 

*(*/) = 1 , 

y{*f) = °- 

The minimum-time problem is subjected to a bound on the control u 

-1 < u < 1 . 



V- AXIS 
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G. Example Problems Solved Using the SGR Technique 



Figure G.1.1 Bounded Brachistochrone Problem 



G.3 A GEODESIC PROBLEM [42; 
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The above problem is the same as the one posed in Chapter (3) of Reference [29 
(pps 112-114). Solution of this optimization problem involves “bang-bang” control, 
i.e. the control “bang” from « =-l to « = +1 and vice versa. 

The inequality constraint on the control can be converted into a path equality 
constraint with the use of an auxiliary control v 

u 2 + t; 2 - 1 = 0. 

Therefore the transformed problem is given by 

y' = tju . 

Once again, since the end-time of the problem is not specified, we have made a 
change in the independent variable like the one used in Example (1). 

The initial conditions of the transformed problem are 

x(0) = 0, 

y(o) = 0, 

and the terminal conditions are : 

x(l) = 1, 

. y(i) = 0. 

The path equality constraint of the transformed problem is 

u 2 + t; 2 — 1 = 0. 

The minimum-time of this problem is found numerically to be 2.0041. This mini- 
mum time compares well with that obtained analytically (cf. [29], tj = 2.0000). 
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G. Example Problems Solved Using the SGR Technique 


The “Bang-bang” control of the problem and the optimal trajectory in the distance- 
velocity phase-plane is given in Figure (G.2.2). 

§G.3 A GEODESIC PROBLEM [42] 


min I = tr 

U 

The equations of motion are 

x = u , 
y = v , 
z — w . 

The initial conditions are : 

x(0) = 1 , 
y(0) = 0 , 

*( 0 ) = 0 . 

The terminal condition are : 

- % ■ 

Here u, u, and w are the components of the velocity vector in a 3-D Cartesian space. 
The idea here is to transfer between end points on an ellipsoid, moving at constant 
velocity, in minimum time without leaving the surface of the ellipsoid. 


Thus, the minimum-time problem is subjected to the following state and control 
path equality constraints 
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Figure G.2.2 Bounded Control of j? Plant 
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G. Example Problems Solved Using the SGR Technique 


-r tr -+• U> 2 = 1 , 


x 2 + 2- + = 1 • 


The first constraint states the fact that the magnitude of the velocity is unity. The 
second constraint requires that the travel be made on the surface of an ellipsoid 
with semiaxes a =1, b —2 and c = 3. 


Since the first constraint of the problem is already a path equality constraint, further 
transformation is not necessary. If we take the time derivative of the second con- 
straint equation and substitute into the resultant equation dynamic relations given 
by the equations of motion, we obtain an alternative form of the path equality 
constraint 


yv zw 
xu -I 1 

4 9 


0 


Therefore the tranformed problem is given by 


x' = tju , 

V = t f v. 
z' = tjW . 

Note the change in the independent variable from t to r for the reason explained 
before. 

The initial conditions of the transformed problem are 

*( 0 ) = 1 , 

V(0) = 0, 
x(0) = 0, 
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and the terminal conditions are 


x(l) = 

y( 1) = 

*( 1 ) = 


1 

_ 2 _ 

\/ 3 ’ 

3 

7T 


In view of the need to stay on the surface of the ellipsoid at all times, one of the 
above three initial and terminal conditions is redundant. 

The path equality constraints of the problem are 

w 2 — 1 = 0 , 


2 2 
u t v 


xv 


yv 


zw 

~9 


= 0. 


The minimum-time of this problem is found to be 2.1443. this minimum time com- 
pares well with that obtained using the quasilinearization technique of References 
[49-50; [t f = 2.1439). 

Optimal trajectories in the x — y. x — z and y — z coordinate planes are given in 
Figure (G.3.3). Note that the first two trajectories are part of an ellipse while the 
last one is a straight line. This result is analogous to the shortest path w great-circle r 
travel on the surface of a sphere. 

§G.4 BOUNDED TIME RATE OF CHANGE OF STATE 


min / 

U 


= U 


The equations of motion are given by 


x = u, 


y = u 2 — r 2 — 0.5 . 
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Figure G.S.S A Geodesic Problem 
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The initial conditions are : 

x(0) = 0, 

V(0) = 0. 

The terminal conditions are : 

*(*/) = 1 , 

This minimum-time problem is subjected to a bound on the time rate of change of 
the state y from below 

y > -0.5 

The inequality constraint can be converted into a path equality constraint by the 
use of the dynamic equation of y. The result is 

u 2 — x 2 > 0 . 

The inequality can be further transformed into a path equality constraint with the 
use of an auxiliary control v 

u 2 — x 2 - v 2 = 0 . 

Therefore the transformed problem is given by 

x' = t f u, 

y' = tf(u 2 — x 2 — 0.5) . 

Note the change in the independent variable from t to r for the reason given before. 

The initial conditions of the transformed problem are : 

*(0),= 0, 

V(0) = 0, 
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G. Example Problems Solved Using the SGR Technique 


and the terminal conditions are : 

x(l) = 1, 

/ v 7T 

*(') = — t - 

The path equality constraint of the transformed problem is 

u 2 - x 2 - v 2 = 0 . 

The minimum-time of this problem is found to be 1.8191, which compares well 
with that computed using the quasilinearization technique of Referene [49-50] and 
reported in Reference [48] {t j = 1.8222). 

The optimal time history of the state and the control of the problem are given in 
Figure (G.4.4). Note that in order to satisfy the bound on the time rate of change 
of y, the control u must at all times be larger or equal to x. 


G.4 BOUNDED TIME RATE OF CHANGE OF STATE 
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Figure G.4.4 Bounded Time-rate-of-change of State Control 






